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ABSTRACT 


Direct-Sequence  Spread-Spectrum  (DS-SS)  is  among  the  preferred  modulation 
techniques  for  military  applications.  DS-SS  offers  three  greatly  desired  characteristics.  It 
allows  for  the  development  of  Low  Probability  of  Detection  (LPD)  and  Low  Probability 
of  Intercept  (LPI)  systems  and  has  a  very  good  performance  in  fading  channels.  This  the¬ 
sis  investigates  the  performance  of  the  “Cross-Product  RV  (CRV)  decomposition”  as  the 
basis  of  blind-equalization  algorithms.  The  CRV  is  a  rank-revealing  decomposition  alter¬ 
native  to  the  Eigenvalue  Decomposition  (EVD)  that  can  provide  a  recursively  updated 
estimate  of  the  signal  and  noise  subspace  at  a  reduced  computational  cost.  The  CRV  up¬ 
dating  algorithm  is  implemented  in  MATLAB  and  evaluated  in  a  previously  proposed 
communication  scheme  intended  for  use  in  an  underwater  acoustic  network  called  Sea- 
web.  The  underwater  channel  is  modeled  with  the  Monterey-Miami  Parabolic  Equation 
Model  (MMPE)  for  various  multipath  perturbations.  The  receiver  performance  is  exam¬ 
ined  using  a  Monte  Carlo  simulation.  Bit-error  rates  versus  signal-to-noise  ratio  are  pre¬ 
sented  for  various,  noise  assumptions,  and  receiver  synchronization  assumptions. 
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EXECUTIVE  SUMMARY 


Underwater  wireless  networks  have  many  military  and  eommereial  applieations. 
With  sueh  networks,  the  eontrol  and  eoordination  of  various  assets,  sueh  as  mini¬ 
submarines,  divers,  submerged  instruments  and  Unmanned  Underwater  Vehieles  (UUVs) 
is  possible  without  the  burden  of  eables  or  human  intervention. 

The  underwater  environment  is  possibly  the  most  ehallenging  environment  for 
wireless  eommunieations.  The  propagation  of  sound  is  greatly  affeeted  by  transmission 
loss  (TL),  ambient  noise,  reverberation,  and  the  overall  variability  of  the  ehannel.  These 
limiting  faetors  impose  various  trade-offs  in  the  design  of  a  reliable  eommunieation  sys¬ 
tem.  A  higher  transmission  rate  implies  greater  bandwidth,  higher  TL  and  intersymbol 
interferenee  (ISI)  but  an  almost  negligible  ehannel  variation  over  a  single-bit  interval. 

Direet-sequenee  spread-speetrum  (DS-SS)  is  one  of  the  preferred  eommunieation 
teehniques  in  military  applieations.  Spread-speetrum  modulation  has  a  very  good  per- 
formanee  in  multipath  environments  and  allows  for  the  development  of  Low  Probability 
of  Deteetion  (LPD)  and  Low  Probability  of  Intereept  (LPI)  systems.  Beeause  of  its  de¬ 
sired  eharaeteristies,  DS-SS  is  well  suited  for  shallow-water  eommunieations. 

This  thesis  investigates  the  performanee  of  DS-SS  with  blind-equalization  in 
simulated  oeean  ehannels.  The  equalization  algorithms  examined  are  based  on  the 
“Cross-Produet  RV  (CRV),”  a  relatively  new  subspaee  deeomposition  algorithm  pro¬ 
posed  as  an  alternative  to  the  Eigenvalue  deeomposition  (EVD).  The  CRV  updating  algo¬ 
rithm  provides  a  reeursively  updated  estimate  of  the  signal  and  noise  subspaee  at  less 
eomputational  eost  than  the  EVD.  The  transmitter  and  reeeiver  struetures  are  imple¬ 
mented  in  MATLAB  and  the  ehannel  impulse  response  is  simulated  with  the  Monterey- 
Miami  Parabolie  Equation  model  (MMPE).  The  impulse  responses  extraeted  account  for 
various  channel  perturbations  (interface  roughness,  internal  waves,  etc.)  and  also  for  the 
Doppler  effect  caused  by  the  source  motion. 

The  performance  of  the  proposed  communication  scheme  is  evaluated  using 
Monte  Carlo  simulation  for  various  signal-to-noise  ratios  (SNR)  and  design  parameters 

xvii 


(packet  size,  data  rate  and  ehannel  length).  For  all  the  individual  perturbations  simulated 
with  the  MMPE  model  a  Bit  Error  Rate  (BER)  below  10'^  is  aehievable  for  a  relatively 
low  SNR.  In  the  ease  where  all  perturbations  are  eombined  or  synehronization  between 
the  reeeiver  and  the  transmitter  is  not  aehieved,  the  reeeiver’s  performanee  is  limited. 


I.  INTRODUCTION 


A.  BACKGROUND 

The  primary  purpose  of  this  thesis  is  to  improve  upon  two  previously  proposed 
Direet-Sequenee  Spread-Speetrum  (DS-SS)  communieation  schemes  [1,2]  intended  for 
use  in  an  underwater  acoustic  communication  network  known  as  Seaweb.  Seaweb  is  be¬ 
ing  developed  by  the  Space  and  Naval  Warfare  Systems  Center,  San  Diego,  [3]  and  is 
motivated  by  a  requirement  for  wide-area  undersea  surveillance  in  littoral  waters.  This 
undersea  wireless  network  provides  the  command,  control  and  communications  infra¬ 
structure  [4]  for  coordinating  appropriate  assets,  such  as  mini-submarines,  divers  and 
Unmanned  Underwater  Vehicles  (UUVs)  to  accomplish  a  given  mission  in  an  arbitrary 
ocean  environment.  The  first  underwater  communication  scheme  [1]  involved  a  Direct- 
Sequence  Differential  Binary  Phase-Shift  Keying  with  quadrature  spreading  (DS-IQ- 
DBSK).  Error-correction  coding  and  a  RAKE  receiver  for  multipath  reception  were  used 
to  improve  system  performance.  The  second  and  latest  communication  scheme  as  part  of 
a  previous  thesis  research  [2]  tried  to  develop  an  improved  receiver  structure  based  on 
blind-equalization  algorithms.  In  both  schemes  [1,  2]  the  impulse  response  of  the  under¬ 
water  acoustic  channel  was  generated  using  the  Bellhop  acoustic  propagation  model  [5]. 
The  Bellhop  model  is  a  “variant  of  Gaussian  beam  tracing  and  is  especially  suited  for 
high-frequency  deep-water  problems  where  normal  mode  and  parabolic  equation  (PE) 
models  are  not  practical  [6]”.  In  this  thesis,  the  channel  impulse  response  is  estimated  us¬ 
ing  the  Monterey-Miami  Parabolic  Equation  (MMPE)  model  described  in  [7].  This  is  a 
range  dependent  wave-theory  model  based  upon  the  parabolic  approximation  of  the  wave 
equation  [6].  It  is  a  widely  used  numerical  tool  for  solving  the  acoustic  wave  equation 
especially  for  low  frequencies  and/or  range-dependent  environments  where  it  is  computa¬ 
tionally  more  efficient. 

B,  SEAWEB  UNDERWATER  NETWORK  REQUIREMENTS 

Eike  all  other  military  applications,  Seaweb  must  ensure  communication  security 
in  two  ways.  It  must  limit  the  ability  of  unauthorized  users  to  detect  the  presence  of  the 
communication  signal  (i.e.,  low  probability  of  detection)  and  at  the  same  time  decrease 
the  ability  of  a  hostile  observer  to  “listen  in”  when  communications  are  taking  place  (i.e.. 
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low  probability  of  intercept).  The  aforementioned  requirements  are  met  by  using  DS-SS 
as  the  modulation  scheme  for  packet  transmission.  DS  spread-spectrum  systems  are  able 
to  hide  the  signal  in  the  background  noise  by  operating  at  low  bit-energy  per  noise- 
spectral-density  Also  due  to  the  pseudo-random  properties  of  the  code  (chip¬ 

ping  sequence)  used  to  spread  the  message  signal,  DS  systems  have  a  very  low  probabil¬ 
ity  of  intercept.  The  message  can  be  demodulated  only  by  the  intended  receivers  who 
know  the  exact  replica  of  the  code  used  at  the  transmitter  [8]. 

Other  Seaweb  requirements  include  utility  packets  of  fixed  length  (72  bits)  and  a 
limited  operating  bandwidth  of  5  kHz  (9-14  kHz  operating  frequencies).  The  desired  in¬ 
formation  transmission  rate  is  40  to  100  bits  per  second  (bps)  at  a  bit-error  rate  (BER)  of 
10'^.  The  communication  range  is  between  three  to  five  kilometers  and  the  operating 
depth  is  on  the  order  of  50  to  200  m. 

C.  GOALS  AND  METHODOLOGY 

The  first  step  in  improving  the  previously  proposed  signaling  schemes  is  to  use  a 
more  accurate  underwater  acoustic  model  to  account  for  the  temporal  and  spatial  variabil¬ 
ity  of  the  channel.  As  already  mentioned,  this  model  is  the  MMPE  model  and  will  be  de¬ 
scribed  in  more  detail  in  Chapter  III  of  this  thesis.  Eor  the  sake  of  comparison  various 
impulse  responses  are  extracted  corresponding  to  different  undersea  environments. 

The  combination  of  the  DS-SS  transmitting  scheme  and  the  receiver  structure 
based  on  blind  equalization  algorithms  was  proven  very  robust  as  part  of  a  previous  thesis 
[2].  However,  the  algorithm  worked  completely  offline  and  the  whole  received  packet 
was  needed  in  order  to  calculate  the  filter  coefficients  and  to  demodulate  the  transmitted 
signal  efficiently.  As  part  of  this  thesis  research  a  relatively  new  subspace  tracking  algo¬ 
rithm,  the  cross  product  RV,  also  known  as  CRV  [9],  is  used  to  allow  for  a  real-time  es¬ 
timation  of  the  desired  subspaces.  The  transceiver  with  the  CRV  algorithm  is  imple¬ 
mented  in  MATEAB  for  various  combinations  of  multipath,  white  and  colored  noise. 

D,  BENEFITS  OF  STUDY 

Underwater  wireless  communications  are  much  desired  in  military  applications. 
The  possibility  to  maintain  signal  transmission  without  the  burden  of  cables  enables  the 
operation  of  underwater  robots  and  vehicles  (UUVs)  and  the  gathering  of  data  from  sub- 
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merged  instruments  without  human  intervention.  An  underwater  mobile  network  such  as 
Seaweb  has  many  military  and  commercial  applications  (submarine  communications, 
ocean  exploration,  etc.). 

E.  THESIS  ORGANIZATION 

The  remainder  of  this  thesis  is  organized  into  five  chapters.  Chapter  II  discusses 
the  general  characteristics  of  underwater  acoustic  communication  channels.  Chapter  III 
provides  a  description  of  the  Monterey-Miami  Parabolic  Equation  (MMPE)  model  and  an 
outline  of  the  theory  behind  the  simulated  perturbations.  Chapter  IV  describes  the  CRV 
algorithm,  proposed  as  a  substitute  of  the  Eigenvalue  Decomposition  (EVD)  used  in  [2]. 
The  CRV  updating  algorithm  can  provide  a  recursively  updated  estimate  of  the  signal  and 
noise  subspace.  Chapter  V  reviews  the  previously  proposed  signaling  scheme  [2]  and 
presents  the  new  scheme  based  on  the  CRV  algorithm.  The  performance  of  the  new 
transmitter/receiver  structure  is  evaluated  under  different  channel  conditions  and  design 
parameters.  Einally,  Chapter  VI  reviews  and  summarizes  the  important  results  and  rec¬ 
ommends  follow-on  work. 
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II.  CHARACTERISTICS  OF  UNDERWATER  ACOUSTIC 
COMMUNICATION  CHANNELS 


The  sea  forms  a  remarkably  complex  waveguide  medium  for  the  propagation  of 
sound.  The  purpose  of  this  chapter  is  to  describe  the  dominant  factors  affecting  sound 
propagation  in  shallow  water  and  the  difficulties  those  factors  impose  on  designing  a  reli¬ 
able  communication  system.  The  observation  of  sound  waves  is  greatly  affected  by 
transmission  loss  (TL),  ambient  noise,  and  the  temporal  and  spatial  variability  of  the 
channel.  Transmission  loss  and  noise  determine  the  signal-to-noise  ratio  (SNR),  the 
maximum  transmission  range,  and  the  available  communication  bandwidth.  The  temporal 
and  spatial  variability  of  the  channel  induced  intersymbol  interference  (ISI),  which  influ¬ 
ences  system  design  (i.e.,  the  choice  of  the  modulation/signal  detection  method  and  the 
transmitter/receiver  structure). 

A,  TRANSMISSION  LOSS  (TL) 

The  weakening  of  an  acoustic  signal  as  it  travels  through  the  sea  is  quantitatively 
described  by  transmission  loss,  TL.  According  to  [10]  and  [1 1],  TL  is  defined  as  the  ratio 
in  decibels  between  the  acoustic  intensity  I{r)  in  W/m  at  a  field  point  (of  distance  r) 

and  the  intensity  measured  at  a  distance  Tq  =  1  m  from  the  sound  source  (transmit¬ 

ter). 


TL  [dB  ref  r^]  =  -10  log 


M 

L 


=  10  log 


I{r) 


(2.1) 


For  both  plane  and  spherical  waves,  the  acoustic  intensity  is  proportional  to  the 
square  of  the  pressure  amplitude;  therefore,  TL  can  also  be  expressed  as 


TL  [dB  ref  ]  =  -20  log 


P{r) 


=  20  log 


P{r) 


(2.2) 


Transmission  loss  can  be  decomposed  into  two  terms,  ,  caused  by  en¬ 
ergy  spreading  and  ,  due  to  sound  attenuation.  Energy  spreading  depends  on 

the  propagation  distance  and  the  underwater  environment  characteristics.  For  the  case  of 
an  omni-directional  sound  source  in  an  unbounded  medium,  sound  energy  undergoes 
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spherical  spreading  and  TL  due  to  spreading  is  proportional  to  the  square  of  the  transmis¬ 
sion  range, 

TLu^bounded, Isospeed  [dB  rcf  1  m]  =  20 log r.  (2.3) 

In  isospeed  underwater  channels  where  the  medium  is  bounded  by  two  plane  and 
parallel  surfaces,  the  TL  in  the  far-field  is  proportional  to  the  first  power  of  the  propaga¬ 
tion  distance  (cylindrical  spreading). 


Bounded,  Isospeed  [dB  rcf  1  m]  =  1  0  lOg  T.  (2.4) 

In  shallow  water  channels  acoustic  energy  undergoes  both  types  of  spreading 
mentioned  above-spherical  spreading  close  to  the  source  and  cylindrical  spreading  at 
great  distances  [12].  Thus  a  more  practical  expression  for  the  TL  in  shallow  water  is 
given  by  [10] 


Shallow  Water  Spreading  [dB  rcf.  1  m] 


20  log  r  r<r^ 

lOlogr +  101ogr,  r  > 


(2.5) 


The  transition  range  (i.e.,  the  range  at  which  the  change  from  spherical  to  cylindrical 

spreading  occurs)  for  the  case  of  an  isospeed  shallow  water  channel  with  fast  bottom,  is 
estimated  by  [10] 


n  = 


(2.6) 


where  H  is  the  bottom  depth  in  meters  and  6^  is  the  critical  angle  in  degrees. 


Figure  1 .  Spreading  of  Acoustic  Energy  in  a  Shallow  Water  Channel. 
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Attenuation  of  the  field  is  the  effeet  of  sound  absorption.  Absorption  involves  the 
proeess  of  eonverting  aeoustie  energy  into  heat,  thus  resulting  in  a  true  loss  of  energy. 
The  dominant  eauses  of  absorption  in  seawater  are  the  ionie  relaxation  of  magnesium  sul¬ 
fate  (MgS04),  the  borie  aeid  (H3BO3)  ionization,  and  the  shear  viseosity  of  seawater  [10, 
11].  Sound  absorption  is  quantitatively  deseribed  by  the  absorption  eoeffieient  a  with 
units  of  dB/m.  One  of  the  most  aeeurate  expressions  for  the  frequeney  dependent  absorp¬ 
tion  eoeffieient  is  provided  by  Fisher  and  Simmons  [10]  and  is  based  on  experimental 
data; 

(a  B  \ 

a[dB/m]=  +  +  C  f\  (2.7) 

U  +/  /2  +/  ) 

The  terms  A,  B,  C,  depend  on  temperature  and  hydrostatie  pressure  whereas  the  relaxa¬ 
tion  frequeneies  y]  and  /j  are  funetions  of  temperature  only.  Speeifioally,  after  the  ap¬ 
proximation  given  in  [10] 


/  [Hz]=780/'"^ 

(2.8) 

/j  [Hz]  =  42000/"'', 

(2.9) 

Z  =  8.3  X 1 0^'  (5  /  35)  ^ 

(2.10) 

(2.11) 

(2.12) 

where  Z  is  the  water  depth  in  km,  T  is  the  temperature  in  °C,  S  is  the  salinity  in  ppt  and 
pH  is  the  aeidity.  The  absorption  eoeffieient  in  seawater  and  frequeneies  between  100  Hz 
and  1  MHz  is  shown  in  Figure  2  (from  Equation  (2.7)  with  parameter  values  given  in  the 
figure  eaption).  In  shallow  water  (Z  <  150  m),  depth,  temperature  and  salinity  have  little 
effeet  on  the  absorption  of  sound. 
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Figure  2.  Absorption  Coefficient  (dB/km)  in  Seawater  (pH  =  8,5  =  35  ppt,  Z  =  0. 1 

km,  T=  15°C). 

Now  that  simple  forms  for  energy  spreading  and  sound  absorption  are  defined,  a 
general  approximation  of  the  transmission  loss  in  shallow  water  is  given  by  [9,  11] 


TL  [dB  ref  Im]  = 


20  log  r  +  ar  f-f~t 

lOlogr +  101ogrj +ar  r  >  r/ 


(2.13) 


where  the  transmission  range  r  and  the  transition  range  r,  are  in  meters  and  the  absorp¬ 
tion  coefficient  a  in  dB/m.  The  range  and  frequency  dependency  of  TL  for  a  fast  bottom 
environment  is  illustrated  in  Figure  3  (for  the  parameters  given  in  the  captions  of  Figures 
2  and  3).  The  frequencies  chosen  to  evaluate  Equation  (2.13)  are  the  lowest  and  highest 
operating  frequencies  of  Seaweb,/;  =  9  kHz  and  f2  =  \4  kHz,  respectively.  A  third  fre¬ 
quency  (/)  =  50  kHz)  much  higher  than  f2  is  used  to  demonstrate  the  logarithmic  increase 
of  the  TL  at  high  frequencies.  In  Figure  4,  the  TL  is  plotted  as  a  function  of  frequency  for 
the  maximum  transmission  range  (5  km)  according  to  Seaweb  specifications. 
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Figure  3.  Range  and  Frequency  Dependency  of  XL  =150m,  6^  =30° 


Figure  4.  XL  for  Frequencies  5  kHz  -  100  kHz  =150m,  6^  =30° 
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As  seen  in  Figure  4,  for  a  frequeney  of  12  kHz  (the  mean  frequeney  of  Seaweb), 
the  XL  is  approximately  45  dB;  whereas  for  a  frequeney  of  50  kHz  the  XL  inereases  loga- 
rithmieally  assuming  the  value  of  110  dB.  Xhe  inerease  of  the  XL  at  high  frequeneies  is 
the  reason  the  available  eommunieation  bandwidth  of  Seaweb  is  limited  to  5  kHz.  Un¬ 
derwater  eommunieation  ehannels  are  severely  band  limited  eompared  to  the  ehannels 
used  in  RF  eommunieations.  Beeause  of  the  limited  bandwidth,  it  is  very  diffieult  to 
aehieve  a  reliable  link  at  high  bit  rates.  By  inereasing  the  data  rate,  the  period  of  the  mes¬ 
sage  bit  deereases  and  the  bandwidth  of  the  signal  inereases.  For  the  sake  of  eomparison, 
eaeh  ehannel  of  the  802.1  Ig  wireless  network  aeeording  to  the  IEEE  standard  has  a 
bandwidth  of  25  MHz  and  a  maximum  throughput  of  54  Mbps,  whereas  an  underwater 
network  ean  only  have  a  few  kHz  of  bandwidth  and  a  mueh  lower  data  rate. 

B,  AMBIENT  NOISE 

Xhe  underwater  environment  everywhere  eontains  ambient  noise.  Ambient  noise 
is  the  total  noise  background  minus  the  noise  generated  by  the  transmitting  equipment 
used  in  a  certain  application.  Noise  observed  in  the  ocean  exhibits  strong  frequency  de¬ 
pendence  as  well  as  geographical  dependence.  In  deep  water,  there  are  many  sources  of 
ambient  noise  and  each  source  dominates  in  a  different  frequency  band.  At  low  frequen¬ 
cies  (below  20  Hz)  the  main  noise  sources  are  the  ocean  turbulence  and  the  various  tec¬ 
tonic  processes  such  as  earthquakes  and  volcanic  eruptions  [13].  In  the  frequency  band  50 
to  500  Hz,  distant  ship  traffic  and  distant  storms  are  the  main  sources  of  noise  [12].  Be¬ 
tween  0.5  and  50  kHz,  the  detected  noise  is  associated  with  the  breaking  of  the  sea  sur¬ 
face  and  bubble  formations.  Xhe  state  of  the  local  sea  surface  is  weather  dependent  and 
closely  connected  with  the  prevailing  wind  speed.  Xhe  greater  the  wind  speed  the  greater 
the  ambient  noise  contribution.  Above  50  kHz,  thermal  noise  due  to  the  molecular  agita¬ 
tion  of  the  water  molecules  predominates  [10].  Xhe  interested  reader  can  read  more  about 
ambient  noise  and  its  sources  in  References  [14]  and  [15].  Xhe  ambient  noise  spectrum 
level  (ANL)  due  to  the  combined  effects  of  the  aforementioned  sources  is  illustrated  in 
Eigure  5.  In  the  frequency  band  (9  to  14  kHz),  at  which  Seaweb  operates,  the  wind  speed 
(or  sea  state)  is  the  main  source  of  noise.  Eor  example,  at  12  kHz  and  a  wind  speed  of  15 
knots,  the  ambient  noise  spectrum  level  is  approximately  45  dB  ref  1  pPa.  In  general,  the 
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ANL  decreases  when  the  frequency  increases,  and  falls  at  a  rate  of  15  dB/decade  when 
the  frequency  range  is  between  0.5  and  50  kHz. 


Figure  5.  Average  Deep-Water  Ambient  Noise  Spectrum  Level  [From  Ref.  12]. 


Wenz  curves  [15]  (or  modified  ones  like  Figure  5)  are  commonly  used  to  estimate 
the  ambient  noise  spectrum  level  under  various  conditions  (e.g.,  shipping  density,  sea 
state,  etc.).  The  following  expressions  [16]  are  a  quantitative  description  of  Wenz  curves: 

ANL„„^,=17-301og(/),  (2.14) 

=  40  +  20(D  -  0.5)  +  26  log(/)  -  60  log(/  +  0.03),  (2.15) 

ANL^.„,  =  50  +  7.5  (w°" )  +  20  log(/)  -  40  log(/  +  0.4),  (2.16) 

ANL„=-15  +  201og(/),  (2.17) 


and 


ANL 


Deep  Water 


=  101og 


ANLr. 


10 


+  10 


^^^Shipping 

io 


ANL^ 


+  10 


10 


+  10 


10 


(2.18) 


V  y 

where  D  is  the  shipping  density  with  values  between  0  and  1  (1  for  heavy  traffic),  w  is  the 
wind  speed  in  m/s,/  is  the  frequency  in  kHz,  and  ANL^^^p  is  the  ambient  noise  level 
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in  deep  water  with  units  of  dB  ref  1  pPa.  Equation  (2.18)  is  evaluated  for  various  ship¬ 
ping  densities  in  Figure  6  and  for  various  wind  speeds  in  Figure  7. 


Figure  6.  Ambient  Noise  Fevel  in  Deep-water  for  Wind  Speed  w=  \5  m/s. 


Figure  7.  Ambient  Noise  Fevel  in  Deep-water  for  Moderate  Shipping  {D  =  0.5). 
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In  shallow  water,  ambient  noise  levels  are  very  difficult  to  define  because  two  ad¬ 
ditional  noise  sources,  man-made  noise  and  biological  noise  (dolphins,  snapping  shrimp, 
etc.)  exist.  These  sources  vary  greatly  both  in  time  and  location.  Man-made  noise  is 
higher  close  to  harbors  and  coastal  areas  with  industrial  activity.  Noise  from  biologies  is 
in  general  unpredictable.  It  is  quite  difficult  to  estimate  beforehand  and  existing  meas¬ 
urements  (if  any)  are  only  valid  in  the  specific  location  they  were  taken.  However,  in  the 
absence  of  industrial  and  biological  noise,  the  ambient  noise  level  depends  mostly  on  the 
strength  of  locally  generated  wind  and  is  about  5  dB  higher  than  the  deep  water  levels  for 
the  same  frequency  [15]. 

Both  transmission  loss  and  the  ambient  ocean  noise  contribute  to  the  degradation 
of  the  acoustic  signal  as  it  propagates  away  from  the  source.  By  combining  Equations 
(2.13)  and  (2.18)  and  assuming  that  the  noise  generated  by  Seaweb  equipment  is  zero,  we 
can  calculate  the  signal-to-noise  ratio  (SNR)  of  the  underwater  communication  channel 
from 

SNR[dB]  =  SL  +  DlT+DlR-TL-ANL,  (2.19) 

where  SL  is  the  source  level  in  dB  ref.  I  pPa,  Im,  and  DI.j.,  DIj^  are  the  directivity  indi¬ 
ces  of  the  transmitter  and  the  receiver,  respectively.  The  directivity  index  is  a  quantitative 
measure  of  the  focusing  of  acoustic  energy,  and  for  an  omni-directional  source  or  re¬ 
ceiver,  DI  =  0 .  The  SLof  the  low  frequency  omni-directional  “telesonar”  transducer  AT- 
408  (a  candidate  for  Seaweb)  according  to  the  OEM  specifications  [17]  is  180  dB  ref  I 
pPa,  Im  at  full  power.  The  resulting  SNR  of  the  channel  is  illustrated  in  Eigure  8.  We  can 
see  from  this  figure  that  for  a  frequency /=  14  kHz  and  a  transmission  range  of  r  =  5  km, 
the  SNR  is  approximately  76  dB.  Needless  to  say,  this  is  only  a  rough  and  rather  optimis¬ 
tic  estimate  of  the  SNR  found  in  a  real  underwater  environment.  The  frequency  depend¬ 
ent  SNR  for  various  transmission  ranges  is  shown  in  Eigure  9.  The  reason  the  frequency 
band  of  9  to  14  kHz  is  chosen  for  the  Seaweb  network  is  clear  since  this  band  accounts 
for  the  flat  portion  (SNR  is  almost  constant)  of  the  curves  corresponding  to  r  =  3  and  r  = 
5  km. 
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Figure  8.  SNR  for  the  Seaweb  Operating  Frequeneies  (wind  speed  15  m/s,  D  =  0.5). 


Figure  9.  SNR  for  Different  Transmission  Ranges  (wind  speed  15  m/s,  D  =  0.5). 
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c. 


TEMPORAL  AND  SPATIAL  VARIABILITY  OF  THE  CHANNEL 


The  characteristics  of  the  underwater  communication  channel  vary  both  in  time 
and  location.  The  transmitted  signal  is  subject  to  multipath  propagation  and/or  Doppler 
spreading.  Together,  these  factors  may  cause  a  severe  degradation  in  the  received  signal 
strength. 

1,  Multipath  Propagation 

The  acoustic  wave  as  it  travels  through  the  ocean  rarely  follows  a  single  path.  Due 
to  reflection  at  the  sea  boundaries,  refraction,  and  volume  scattering  within  the  water  col¬ 
umn  or  the  sub-bottom  structure  caused  by  random  inhomogeneities,  multiple  echoes  of 
the  transmitted  signal  arrive  at  the  receiver.  These  time-delayed  echoes  vary  both  in  am¬ 
plitude  and  phase  and  interfere  with  one  another.  This  effect  is  called  multipath  propaga¬ 
tion  or  time  spreading  since  travel  time  is  different  for  each  path.  Multipath  characteris¬ 
tics  are  different  in  shallow  and  deep-water  channels.  In  shallow  water,  the  multipath  is 
severe  because  of  repeated  reflections  from  the  sea  surface  and  bottom.  Each  reflection 
causes  an  attenuation  of  the  sound  wave  prohibiting  long-range  communications  without 
a  significant  loss  of  acoustic  energy  [13].  The  relative  delay  time  between  the  first  and 
/-th  arrival  of  the  transmitted  signal  at  the  receiver  is  called  excess  delay  T^ .  The  maxi¬ 
mum  excess  delay  between  the  first  arrival  and  the  arrival  of  the  last  distorted  rep¬ 
lica  of  the  original  signal  is  much  greater  in  an  acoustic  communication  channel  than  in  a 
land-based  RE  channel.  This  is  mainly  because  sound  in  water  travels  slower  (~  1500 
m/s)  than  the  electromagnetic  waves  travel  in  the  air  (3x10*  m/s).  In  the  frequency  do¬ 
main,  the  bandwidth  over  which  the  channel  passes  all  spectral  components  with  ap¬ 
proximately  equal  gain  is  called  the  coherence  bandwidth  [18]  and  is  a  useful  tool  for 
comparing  different  multipath  channels: 

Bc~—.  (2.20) 

T 

max 

2,  Doppler  Spreading 

The  underwater  environment  is  time  variant.  The  medium  itself  is  constantly 
moving  due  to  surface  and/or  internal  waves  and  the  source  and  receiver  platform  may 
also  move.  In  shallow  water,  surface  scattering  is  anticipated  to  be  the  most  important 


15 


contributor  to  the  overall  time  variability  of  the  channel.  Surface  scattering  is  caused  by 
the  weather  dependent  roughness  of  the  sea  surface.  The  frequency  of  an  acoustic  wave 
incident  upon  the  moving  surface  of  the  ocean  is  affected  by  spectral  broadening.  This 
phenomenon  is  called  Doppler  spread.  The  Doppler  spread  for  the  case  of  a  sound  wave 
of  frequency  /  reflected  at  the  sea  surface  is  given  by  [13] 


B 


D 


(1.75x10^') 


^  f  cos  6*^ 

I  c  J 


(2.21) 


where  is  the  Doppler  spread  in  Hz,  w  is  the  wind  speed  in  m/s,  6*  is  the  grazing 

angle,  and  c  is  the  speed  of  sound  (  ~  1500  m/s).  For  a  transmission  range  much  greater 
than  the  chaimel  depth  (z  »  r  ),  6  is  very  small  and  cos  6*  » 1 .  In  Figure  10,  the  Doppler 
spread  for  a  frequency  /  =  14kHz  and  wind  speed  >v  =  15m/s  (sea  state  3) 
is  approximately  9  Hz  and  the  frequency  spectrum  of  the  reflected  wave 
is/ ±5^  =14 ±0.009  kHz. 
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Figure  10.  Doppler  Spread  Due  to  Surface  Reflection. 
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The  second  parameter  used  to  describe  the  time  variability  of  the  channel  is  the 
coherence  time  ,  which  is  inversely  proportional  to  the  Doppler  spread  (for  simplicity 


the  Doppler  spread  is  considered  to  be  the  only  factor  contributing  to  coherence  time), 

7f=— .  (2.22) 

Bo 

The  meaning  of  coherence  time  is  that  two  signals  arriving  at  the  receiver  with  a  time 
separation  less  than  will  be  affected  by  the  channel  in  the  same  way  [18]. 

The  relation  between  the  aforementioned  parameters  B^,  B^,  that  de¬ 
scribe  the  overall  variability  of  the  channel  leads  to  four  distinct  types  of  fading  that  are 
summarized  in  Figure  1 1 ; 


Figure  11.  Different  Types  of  Fading  Channels  [After  Ref  18]. 


The  first  two  types  of  fading  are  frequency  selective  and  frequency  non-selective 
fading.  A  frequency  non-selective  channel  implies  that  the  coherence  bandwidth  of  the 
channel  B^  is  much  greater  than  the  bandwidth  of  the  transmitted  signal  or  equiva¬ 
lently  the  maximum  excess  delay  is  much  less  than  the  digital  symbol  (bit)  pe¬ 
riod?],  .  In  this  case,  which  is  the  best  for  wireless  communications,  the  channel  induces 
no  intersymbol  interference  (ISI).  On  the  other  hand,  when  frequency  selective  fading 
(Bc<  B^  )  occurs,  then  T,  <  and  ISI  is  unavoidable. 

The  other  two  types  of  fading  are  based  on  the  coherence  time  of  the  channel.  The 
first  one  is  called  a  “fast-fading”  channel  and  the  second  one  a  “slow-fading”  channel.  In 
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a  fast-fading  communication  link  the  bit  duration  is  greater  than  the  eoherence  time  of  the 
channel.  The  ehannel  impulse  response  ehanges  rapidly  during  the  time  it  takes  to  trans¬ 
mit  one  symbol.  In  this  type  of  fading  a  severe  distortion  of  the  signal  oecurs  and  a  reli¬ 
able  communieation  link  cannot  be  maintained.  In  a  slow-fading  ehannel, 

{B,  «  Bg  ),  and  less  distortion  is  introduced  as  a  result  of  the  wave  propagation  in  the 
water. 

In  this  chapter,  we  have  seen  the  reason  the  shallow-water  channel  is  eonsidered 
to  be  one  of  the  most  challenging  ehannels  for  wireless  eommunications.  The  system  en¬ 
gineer  must  take  into  aceount  all  the  limiting  factors  (TL,  ambient  noise,  fading)  imposed 
by  the  channel  in  order  to  design  a  reliable  eommunication  link.  Beeause  of  the  time- 
varying  multipath,  there  is  always  a  trade-off  in  choosing  the  transmission  rate.  A  higher 
bit  rate  implies  greater  bandwidth,  higher  transmission  loss  and  ISI  but  an  almost  negligi¬ 
ble  ehannel  variation  over  a  single -bit  interval.  In  order  to  undo  the  effects  caused  by  the 
channel,  a  statistical  model  of  the  channel  impulse  response  is  needed.  In  the  next  chap¬ 
ter,  we  describe  the  steps  needed  to  extract  the  ehannel  impulse  response  by  using  one  of 
the  existing  ehannel  models  called  MMPE,  which  is  based  on  the  parabolic  approxima¬ 
tion  of  the  wave  equation. 
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III.  MODELLING  OF  THE  UNDERWATER  CHANNEL 


In  order  to  allow  for  the  implementation  of  an  effieient  eommunieation  seheme,  a 
statistieal  representation  of  the  underwater  ehannel  is  needed.  In  underwater  aeousties 
various  numerieal  methods  are  used  to  deseribe  the  propagation  of  sound.  All  eurrent 
methods  solve  the  wave  equation  for  a  speeifie  environment  by  following  two  different 
theoretieal  approaehes.  The  first  approaeh  is  based  on  wave  theory  (normal-mode,  para- 
bolie  equation  models)  and  the  seeond  one  on  ray  aeousties  (ray  models)  [12].  Eaeh  the¬ 
ory  has  its  own  advantages  and  disadvantages.  Wave  theory  models  give  a  eomplete  solu¬ 
tion  and  are  best  suited  for  shallow-water  and  low-propagating  frequeneies.  Ray  models 
provide  a  praetieal  visualization  of  the  propagation  of  sound  in  spaee  and  give  their  best 
results  at  high  frequeneies  [11]. 

A,  MONTEREY-MIAMI  PARABOLIC  EQUATION  (MMPE)  MODEL 

In  this  thesis  the  ehannel  impulse  response  is  predieted  using  the  Monterey-Miami 
parabolie  equation  (MMPE)  model  [7].  The  MMPE  is  a  range-dependent,  wave-  theory 
model  based  upon  the  parabolie  approximation  of  the  wave  equation.  Sinee  the  purpose 
of  this  researeh  is  not  the  model  itself,  only  an  outline  of  its  theoretieal  foundation  is  pre¬ 
sented.  We  begin  by  representing  the  time  harmonie  aeoustie  pressure  field  in  eylindrieal 
eoordinates  assuming  azimuthal  symmetry, 

P{r,z,(i>t)  =  p{r,z)e“"‘.  (3.1) 


Next  we  substitute  Equation  (3.1)  into  the  linear,  homogeneous  wave  equation  in  eylin- 
drieal  eoordinates. 


V^p{r,z)--  ^ 


c(r,  z)  dt 

where  c  [m/s]  is  the  speed  of  sound  in  water  and  the  Eaplaeian  defined  as 


V" 


a"  lA  ^ 

dr^  r  dr  dz^ 


(3.2) 


(3.3) 
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This  leads  to  the  time-independent  Helmholtz  equation, 


V V(^, z)  +  ^  .2  p{r, z)  =  0, 

c(r,zf 

whieh  ean  be  faetored  by  introdueing  the  operator  defined  as 

Qop  =  V/^  +  ^  +  1, 

with 


s  =  n^  -  \, 


1  a" 

kl  dz^  ■ 


(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 


In  the  above  expressions,  n  is  the  index  of  refraetion,  Cq  is  the  referenee  sound 
speed  of  the  oeean  volume  (usually  taken  as  1500  m/s),  and  -  col is  the  referenee 

wavenumber.  Note  that  the  eoordinate  dependeney  of  n  and  c  is  dropped  for  simplieity 
but  is  always  implied. 

If  we  take  into  aeeount  the  eylindrieal  spreading,  whieh  is  typieal  in  shallow  water 
and  properly  faetor  the  Helmholtz  equation,  the  out-going  aeoustie  pressure  may  then  be 
defined  as  [19], 

p(r,z)  =  P^^Q;p\(r,z)e'"°\  (3.9) 

where  Pq  is  the  pressure  amplitude  at  distanee  r  =  Rq,  and  z)  is  the  envelope  fune- 

tion,  or  PE  field  funetion.  After  substituting  Equation  (3.9)  into  (3.4)  the  parabolie  equa¬ 
tion  for  the  PE  field  funetion  is  formulated; 

=  -ik,ip  +  ik^Q^pip  =  -ik^H^^y/,  (3.10) 
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where 


(3.11) 

is  an  operator,  defining  the  evolution  of  the  PE  field  funetion  in  range.  By  integrating 
Equation  (3.10)  and  assuming  that is  eonstant  over  a  distanee  Ar,  the  relationship 

between  the  values  of  y/  at  different  ranges  may  be  defined  as 


\j/{r  +  Ar)  =  e^'^'''»"'"’^''V(r),  (3.12) 

where  e  is  a  propagator  that  aeeounts  for  the  ehange  in  the  value  of  the  PE  field 

funetion  y/  with  range.  This  propagator  is  eomputed  using  a  split-step  Eourier  (PE/SSE) 
method  deseribed  in  [20].  The  disere te  fast  Eourier  transform  (EET)  subroutine  allows  the 
PE/SSE  implementation  to  be  represented  by  [7] 


^(r  +  Ar,  z)  =  e 


(r+Ar,z) 


EET< 


-ik„hrT  (k^ ) 


lEET 


y/{r,z) 


(3.13) 


The  and  operators  are  defined  as 


Ua=-{n-\), 


(3.14) 


and 


T  =1- 

op 


1- 


yKj 


(3.15) 


where  is  the  vertieal  wavenumber.  The  output  of  the  model  is  the  PE  field  funetion  y/ 

from  whieh  other  aeoustieal  quantities  ean  be  eomputed.  Eor  example  by  substituting 
Equation  (3.13)  into  (3.9)  we  obtain  the  oeean  frequeney  response, 


G{f)  =  p{r,z)  =  -^y/{r,  z)e'^»'' . 
yjr 


(3.16) 


Eor  broadband  signals  the  MMPE  model  must  be  run  for  eaeh  frequeney  in  the 
band  in  order  to  eompute  the  frequeney  response  at  a  given  range  and  depth.  In  the  fol¬ 
lowing  seetions  of  this  ehapter,  the  theory  behind  the  various  perturbations  modeled  in 
this  thesis  is  deseribed  and  the  results  of  the  simulations  are  presented. 
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B, 


INTERFACE  ROUGHNESS  PERTURBATION 

In  order  to  define  the  interfaee  roughness,  we  assume  a  2-D  interfaee  spectrum  of 


the  form  [21] 

<\K)  = - - - J.  (3.17) 

(i+0;)= 

where  is  the  horizontal  spatial  wavenumber  vector  given  by 

k^=^K^+L\  (3.18) 

K  and  L  are  the  horizontal  wavenumbers  in  the  x-  and  y-directions,  respectively,  Lcorr  is  a 
correlation  length  scale,  (3  is  the  spectral  exponent  and  //  is  a  normalization  factor  de¬ 
fined  in  terms  of  the  root-mean-square  (rms)  roughness  cr  [22], 

\(  0  , 

^  =  -  ^  ^"orr-  (3.19) 

;r  V  2 

'Pox  the  determination  of  the  interface  roughness  affecting  forward  propagation, 
the  one-dimensional  (1-D)  spectrum  along  the  x-axis  (kJ  is  required.  This  is  com¬ 
puted  by  taking  the  1-D  transform  of  the  2-D  interface  spectrum  =  JV2^\K,L) 

along  a  slice  at  y  =  0,  defined  by 

=  (3.20) 

In  cylindrical  coordinates  Equation  (3.20)  yields 

tFy(K)  =  rc^%^Jl  +  Ll^KyPi  (3.21) 

where  y  is  given  by 


(3.22) 


and  r  is  the  Euler  gamma  function.  In  order  to  illustrate  the  interface  roughness  pertur¬ 
bation  to  the  underwater  environment,  the  TE  vs.  depth  and  range  is  computed  for  the 
frequency /=  1 1.5  kHz,  a  bottom  depth  of  150  m,  and  a  source  depth  of  75  m.  The  results 

are  shown  in  Eigure  12  for  an  rms  roughness  value  of  2  m.  By  looking  at  the  two  plots 
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two  observations  can  be  made.  For  the  case  of  the  unperturbed  environment  (top  plot)  the 
down-refracting  sound  speed  profile  causes  significant  interactions  with  the  seafloor.  The 
upper  figure  also  displays  a  clear  interference  pattern  due  to  the  coherent  summation  of 
modes.  On  the  other  hand,  when  the  interface  roughness  is  added,  the  bottom  interactions 
with  the  rough  boundary  create  a  diffuse  field  due  to  the  incoherent  scatter  of  energy  and 
the  breakdown  of  coherence  in  the  mode  interference  pattern. 


Transmission  Loss  (dB  re  1m)  -  No  Perturbation 


Range  (km) 


Transmission  Loss  (dB  re  1m)  -  interface  Roughness  Perturbation 


Range  (km) 


Figure  12.  Interface  Roughness  Perturbation  Compared  to  the  Unperturbed  Environ¬ 
ment  (Source  Depth  =  75  m,  /=  1 1.5  kHz,  Range  0  to  1  km  ). 
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C.  BOTTOM  VOLUME  PERTURBATION 

1.  Volume  Sound  Speed  Fluctuations 

For  the  computation  of  the  volume  sound  speed  fluctuations  affecting  the  forward 
propagation,  the  two-dimensional  (2-D)  vertical  spectrum  of  the  volume  fluctuations 

W^P(ky)  is  required,  where  ky  =  +M^  is  the  wavenumber  in  the  (x,z)  plane.  This 

can  be  determined  by  first  modeling  the  sediment  volume  sound  speed  perturbations  with 
help  of  the  three-dimensional  (3-D)  volume  spectrum  (A,  T,M)  defined  by  [23] 

wP{K,L,M)  =  ^^^{a^K^  +  (3.23) 

In  ^  ^  ’  ' 

where  M  is  the  vertical  wavenumber,  B  is  the  spectral  strength  constant,  and  A  is  the 
horizontal-to-vertical  aspect  ratio  used  to  describe  the  anisotropy  of  fluctuations  in  the 
sediment. 

The  2-D  vertical  volume  spectrum  W^p{K,M)  in  the  (x,z)  plane  is  computed  by 
integrating  (3.23)  over  L  as  follows  [22]; 

wPiK,M)=]  wPiK,L,M)dL  =  ^^^°j[ApK^+Lp  +  M^J^~^dL.  (3.24) 

-00  ^0 

For  the  specific  case  of  =  2,  A  «  5,  and  5  «  5  x  10  [23],  Equation  (3.24)  yields  [22] 

CO  3 

W}P(K,M)=  ^  wP{K,L,M)dL  =  a'[25K^  (3.25) 

-00 

where  a' =  1.25x10  \ 

2,  Volume  Density  Fluctuations 

The  variability  in  bottom  volume  density  p  is  treated  through  the  effective  index 
of  refraction  n'  defined  by  [7] 

11  3/^1  V 

n'^=n^+^  -y^p--  -yp  ,  (3.26) 

2^0  [P  2^/7  ; 

where  k^  is  the  reference  wavenumber  and  n  is  the  index  of  refraction.  By  considering 

the  fact  that  the  sediment  properties  are  horizontally  stratified,  and  assuming  that  the  en¬ 
vironment  is  range-independent  over  a  range  step  Ar  ,  Equation  (3.26)  simplifies  to 
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,2  2,1  1  av  1  dp 

n  =n  H - ^ - ^ - - 

2kl  p  dz^  2yp  dz 


The  definition  of  the  index  of  refraetion  in  the  bottom  is 


(3.27) 


"b  2  ’ 


(3.28) 


where 


The  faetor  Sc/^  is  given  by 


Ca  +  =  +Sc,. 


^Ch=Ci,^{K+^,), 


(3.29) 


(3.30) 


where  is  the  mean  bottom  sound  speed  at  the  interfaee,  Si  is  the  zero-mean  random 
perturbation,  and  =  gjc,^^  is  the  normalized  gradient  of  the  bottom  sound  speed. 

The  relative  fluetuations  in  bottom  volume  density  are  related  with  the  relative 
fiuetuations  in  sound  speed  through  the  expression  [23] 


where 


^  ^  '^{Pr-Po)  ^ ^  ^  (2y)  — 

Po  ~  Pr  J  ^0 


Co  =C4„(1  +  4), 


(3.31) 


(3.32) 


Pr-Po 

2P0-Pr 


(3.33) 


In  this  ease,  p^  =  2650kg/m^  is  the  density  of  the  grain,  and  p^,  Cq  are  the  average  val¬ 
ues  of  the  density  and  sound  speed  in  the  sediment,  respeetively.  From  Equation  (3.31), 
we  obtain  for  the  bottom-volume  density  with  perturbation. 


p=p,+Sp  =  p,+p,{2r)—  =  pJi+2r—  . 


(3.34) 


The  first  and  the  seeond  partial  derivatives  of  (3.34)  with  respeet  to  depth 
z  (depth  gradients  in  or  are  negleeted)  yield 


25 


(3.35) 


dp  lyp^  d 


dz 


Cg  dz 


and 


d^p  _2yPo  d" 
dz^ 


Cq  dz' 


[dc). 


[Sc). 


(3.36) 


These  are  then  ineorporated  into  the  definition  of  the  effeetive  index  of  refraetion,  Equa¬ 
tion  (3.26),  within  the  bottom  volume.  The  bottom-volume  perturbation  in  terms  of  the 
TL  vs.  depth  and  range  is  displayed  in  Figure  13.  As  shown  in  this  figure,  the  seattered 
acoustie  energy  at  short  distances  from  the  source  is  penetrating  into  the  bottom.  How¬ 
ever,  this  has  only  a  minor  impact  in  the  resulting  interference  pattern  compared  to  the 
one  created  by  the  unperturbed  environment. 


Transmission  Loss  (dB  re  1m)  -  Bottom  Voiume  Perturbation 


Range  (km) 


Figure  13.  Bottom  Volume  Perturbation;  TF  vs.  Range  and  Depth  (Source  Depth  = 
75  m,  /=  11.5  kHz,  Range  0  to  5  km  ). 

D.  TURBULENCE  PERTURBATION 

Another  mechanism  affecting  the  coherence  structure  in  an  underwater  environ¬ 
ment  is  the  small-scale  fluctuations  of  sound  speed  within  the  water  volume  caused  by 
turbulence.  The  turbulent  perturbation  field  is  approximated  as  a  random  realization  of 
perturbations  that  are  based  on  an  accepted  power  spectral  density  of  the  turbulence. 


26 


The  variance  of  the  perturbations  to  the  index  of  refraction  ju{r)  is  given  by  [22] 


1  1 

var(//(r))  =  -^  f”  f”  f*  S^{K,L,M)dKdLdM=-^\  S^{k)d^k,  (3.37) 

(2;r)  J-®-'-*-'-®  ^  (2;r)  ^ 

where  S^{k)  =  S^{K,L,M)  is  the  spectral  density  of  the  perturbations  and  k  =  [K,L,M) 

is  the  3-D  wavenumber  vector.  If  the  variance  is  normalized  according  to 

var(//(r))  =  /|//(r)|^\  =  — ^ - \  u\r)d^r,  (3.38) 

\  /  volume*'^ 

where  (  )  denotes  averaging  ,  then  by  combining  Equations  (3.37)  and  (3.38)  we  obtain 

I  (3.39) 

V  ^ )  —00 

Now,  if  F^{k)  is  the  perturbation  transform,  the  index  of  refraction  perturbation  ju{r)  can 
be  written  as 


juir) 


1 

— ^  f  F(k)e‘^-'^d^k, 


(3.40) 


and  its  spatial  autocorrelation  is  given  by 


00  .*  00  oo 

I  ^  j  ]  I  I  Ffk)e“-'d’lc  ||  |  F^d'je 


-00  I  1-00 


■d^r.  (3.41) 


For  F  =0,  Equations  (3.39)  and  (3.41)  yield  the  desired  relation  between  the  signal  spec¬ 
trum  F^(^)  and  the  3-D  spectral  density  of  the  perturbations  S^(k)  [22]; 


F^(k)  =  volume  X  5*^  (^)J  .  (3.42) 

Similarly,  along  the  2-D  plane  of  propagation  (r  =  (x,0,z)) ,  the  relation  between  the  2-D 
spectral  density  V^{K,M),  and  the  2-D  spectrum  of  the  perturbation  field  G^{K,M), 
can  be  written  as  [22] 


G^{K,M)\  =  [Area  xV^{K,M) 


(3.43) 


The  turbulent  perturbation  field  is  approximated  through  the  use  of  the  method  of 
structure  functions  [24],  in  which  large-scale  perturbations  are  not  accounted  for.  The 
three-dimensional  spectral  density  function  E{k),  can  be  expressed  as  [25] 
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2/3,  -11/3 


(3.44) 


E{k)  =  As^^^k- 

where  +l}  ,  s  is  the  energy  dissipation  rate,  and  is  a  scalar  mul¬ 

tiplier.  In  order  to  prevent  E{k)  from  growing  unbounded  when  the  value  of  k  is  small, 
a  lower  wavenumber  threshold  k,  is  defined  so  that  the  A:  dependence  is  [22] 

1 


E{k): 


(k^+K) 


2\ll/6  • 


(3.45) 


Also,  by  applying  the  high-pass  filter 


HPF  = 


(3.46) 


k^+ky 

the  3-D  spectral  density  function  E(k)  is  forced  to  assume  the  value  of  zero  when  A:  =  0  . 
Subsequently,  a  high-frequency  cutoff  is  established  expressed  as  [25] 


Rj^{k)  =  exp 


-q 


yksj 


(3.47) 


where  R^k)  is  the  Batchelor  spectrum,  q  is  an  order  unity  factor,  and  kg  is  the  Batchelor 


wavenumber  defined  by  [22] 


kg  = 


f  _  ^1/4 
8 


^VKtJ 


(3.48) 


In  the  above  expression,  s  [W/kg]  is  the  depth  averaged  kinetic  energy  dissipation  rate. 
Kg  «  lO^^m^/s  is  the  thermal  diffusivity  of  the  sea  water,  and  v  =  1.4x10 m^/s  is  the 
kinematic  viscosity.  In  this  theoretical  development,  s  is  given  by  [26] 


_  0.3 

8  = - 

// 


2 

cosh  ' 

r  N{zy\ 

\L^  ' 

K  fl  jj 

(3.49) 


where  f,  [rad/s]  is  the  inertial  frequency,  «  0.5  m  is  a  measure  of  internal  wave 
intensity,  y,  is  the  mode  number,  and  N{z)  [rad/s]  is  the  buoyancy  frequency  expressed 
in  terms  of  density  p{z)  and  potential  density  Pp{z), 

1 


N(z)  =  g 


1/2 


p(z)  dz 


(3.50) 
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By  following  the  theoretical  approach  described  above,  the  three-dimensional  spectral 
density  S^{k;z)  can  now  be  expressed  by  [22] 


S,ik;z) 


{k^+kff 


R,{k), 


(3.51) 


where  is  the  turbulent  strength  parameter.  In  this  implementation,  ^4^  will  be 

adjusted  manually  to  produce  a  specific  desired  rms  fluctuation. 


The  desired  2-D  spectral  density  is  now  formulated  by  integrating  over  T  (y  =  0) 


only  the  portion  of  the  3-D  spectral  density  that  varies  with  the  turbulent  spectrum.  This 
portion  is  given  by  [22] 


Thus,  the  2-D  spectral  density  (h2^( A:)  can  be  derived  from 


1  1  1  ^ 

{K,M)  =  —  \  ^^dL  =  —  f 


In  +L^ +M^ +ky^^'' 


^  00  ^  ^  GO 

^2^  I 


1  ^ 
a 


11/6 


where  +M^  +  kf  .  With  help  of  integral  tables.  Equation  (3.53)  yields 


0,^{K,M)  = 


1  ^ 


2na^'^  m) 


(3.52) 


(3.53) 


(3.54) 


By  reapplying  the  terms  that  do  not  vary  with  the  turbulent  spectrum,  the  resulting  2-D 
spectral  density  02jj{k)  is  defined  by  [22] 


Kik;z)=  ^ 


r(l)1  ^ 


2n)[T{^))(k^+kfr[k^+k!  I  ^ 


R,{k). 


(3.55) 


The  perturbation  to  the  underwater  environment  because  of  turbulence  is  illustrated  in 
Figure  14  for  an  rms  value  of  2  m/s.  The  breakdown  in  the  coherent  structure  of  the 
modes  is  obvious,  resulting  in  a  highly  diffuse  field  with  random  patterns  of  ray-like 
propagation  (yellow,  high-energy  paths  in  Figure  14). 
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Transmission  Loss  (dB  re  1m)  -  Turbulence  Perturbation 


Range  (km) 

Figure  14.  Turbulence  Perturbation;  TL  vs.  Range  and  Depth  (Source  Depth  =  75  m, 

/=  11.5  kHz,  Range  0  to  5  km). 


E.  INTERNAL  WAVE  PERTURBATION 

Larger-scale  fluctuations  in  sound  speed  may  also  be  generated  from  internal 
waves.  The  internal  wave  perturbation  accounts  for  the  vertical  movement  of  water 
masses  inside  the  sea  body  due  to  buoyancy  forces.  In  this  thesis  a  combination  of  a  lin¬ 
ear  (sinusoidal)  and  a  nonlinear  soliton  internal  wave  perturbation  to  the  sound  speed  is 
modeled.  The  linear  internal  wave  perturbation  is  expressed  as  a  sum  of  sinusoids  [27], 


Ac,i„(r,z)  =  C--exp 


B 


IZcos(Ar), 


(3.56) 


(=1 


where  ACj,i„  (r,z)  is  the  sound  speed  perturbation  at  a  given  range  r  and  depth  z 


B  =  15  m,  A.  =  Inj X. ,  1.  =  1500-300(z -l)  m,  and  C  is  defined  in  such  a  way  that 
the  maximum  value  of  Ac^i^(r,z)  is  1.5  m/s. 


The  nonlinear  soliton  internal  wave  perturbation  Ac^^,!  {r,z^  is  defined  as  [27] 


Ac,„,(r,z)  =  c|exp 


B 


B 


!=1 


sech 


R,-r 

V  A  y 


(3.57) 


where  B  =  25  m,  A.  =  10exp(-0.3(z -l)],  7?j  =  4800m,  7?.^  =  7?.^  -200(7-z)  m, 

D,  =  73430071”  m,  and  C  is  defined  such  that  the  maximum  value  of  the  sound  speed 


30 


perturbation  is  12.5  m/s .  The  internal  wave  perturbation  in  terms  of  the  TL  vs.  depth  and 
range  is  shown  in  Figure  15  (range  0  to  1  km)  and  Figure  16  (range  4  to  5  km).  At  short 
ranges  the  perturbation  is  mostly  due  to  the  linear  internal  waves  and  the  modes  seem  to 
remain  eoherent.  At  long  ranges  from  the  souree  (Figure  16)  the  nonlinear  soliton  internal 
waves  eause  a  much  stronger  perturbation  to  the  field  resulting  in  a  more  random  like 
dispersion  of  the  energy. 


Transmission  Loss  (dB  re  1m)  -  No  Perturbation 


Range  (km) 


Transmission  Loss  (dB  re  1m)  -  internai  Wave  Perturbation 


Range  (km) 

Figure  15.  Internal  Wave  Perturbation  Compared  to  the  Unperturbed  environment 
(Source  Depth  =  75  m,  /=  1 1.5  kFIz,  Range  0  to  1  km). 
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Transmission  Loss  (dB  re  1m)  -  No  Perturbation 


I 

I 


Ran^c  (km) 


lal  Wave  Perturbation 


Figure  16.  Internal  Wave  Perturbation  Compared  to  the  Unperturbed  Environment 
(Source  Depth  =  75  m,  /=  1 1.5  kHz,  Range  4  to  5  km). 

F.  DOPPLER  PERTURBATION 

The  influence  of  Doppler  due  to  source/receiver  motion  is  of  great  interest  in  un¬ 
derwater  acoustic  communications.  Because  of  the  Doppler  effect  the  frequency  of  a 
transmitted  signal  suffers  from  spectral  broadening.  The  actual  frequency  (6*) ,  propa¬ 
gated  from  a  source  that  is  moving  with  speed  ,  can  be  expressed  as  [28] 


/r 


/c<s:l 


1- 


fu  ^ 

I  c  j 


cos 


(^-4) 


fr 


1  + 


fu  ^ 
V  c  J 


cos 


(^-4) 


(3.58) 
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where  is  the  transmitted  frequeney,  c  is  the  sound  speed,  6  is  the  angle  of  propaga¬ 
tion  (relative  to  the  horizontal),  and  is  the  angle  of  souree  motion.  As  the  sound  trav¬ 
els  through  the  medium  it  follows  multiple  paths.  Therefore,  even  a  single  transmitted 
frequeney  will  result  in  multiple  reeeived  frequeneies  at  the  reeeiver. 

If  the  bandwidth  BWj.  of  the  transmitted  signal  is 

BWr=fr  -fr  ■  ,  (3.59) 

T  JT.max  J  T  ,mm.  ^ 

where  fj  and  are  the  minimum  and  maximum  frequeney  eomponents  of  the 

transmission,  then  the  Doppler  shifted  bandwidth  for  forward  propagation  is  given  by 
[28] 

BW"  =BW,  +/r_  |sm((«,)|),  (3.60) 

or 

BWr  =5»'r+^(/,,„.|sin((«,)|  +  A_).  (3.61) 

Equation  (3.60)  accounts  for  a  source  motion  toward  the  receiver  {r  +  )  whereas  (3.61)  is 
for  a  source  motion  away  from  the  receiver  (r  -  ). 

If  the  receiver  is  moving  also,  a  similar  theoretical  approach  yields  [28] 

BW"  +/„.  |sm(«(,)|),  (3.62) 

or  for  the  case  when  the  receiver  is  moving  toward  the  source 

BW'-  =S»;+^^(/,„|sin(A)|  +  /,„.),  (3.63) 

C  ^  / 

where  is  the  receiver  speed,  is  the  angle  of  receiver  motion,  fp  and  fp  are 
the  minimum  and  maximum  frequency  of  propagation,  and  BWp  is  the  bandwidth  of 
propagation  given  by 

=/..„ax (3.64) 
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G.  RESULTS 

The  results  of  the  various  perturbations  deseribed  above  are  presented  in  this  see- 
tion.  The  speeifie  parameters  defining  the  eharaeteristies  of  the  underwater  environment, 
the  souree/receiver,  and  the  resolution  of  the  ealculations  are  summarized  in  Table  1.  For 
each  simulation  the  MMPE  model  was  run  over  a  2048-Hz  bandwidth  (1024  frequen¬ 
cies),  and  the  channel  impulse  response  was  calculated  for  a  source-receiver  distance  of  5 
km  and  a  source  depth  of  75  m  as  in  a  typical  Seaweb  setup.  For  all  modeled  perturba¬ 
tions,  the  transmission  loss  is  presented  as  a  function  of  depth  and  frequency  (Figure  17 
through  Figure  23  -  top  plot).  The  magnitude  of  the  channel  impulse  response  at  a  depth 
of  75  m  versus  the  travel  time  of  the  acoustic  wave  is  illustrated  in  Figure  17  through 
Figure  23  (bottom  plot).  The  impulse  responses  are  normalized  by  the  magnitude  of  the 
unperturbed  response  such  that  its  values  lie  between  -land  1. 

From  the  results  presented  in  the  following  figures  some  interesting  observations 
can  be  made.  The  typical  excess  delay  between  the  first  and  last  arrival  of  the  trans¬ 
mitted  signal  is  about  0.15  seconds  (150  ms).  For  a  symbol  duration  =  0.42ms  corre¬ 
sponding  to  a  chipping  rate  of  2400  chips-per-second  (used  in  the  DS-SS  implementation 
described  in  Chapter  V)  channel-induced  ISI  will  occur  (7^  <  ). 

Among  the  individual  perturbations  examined,  the  one  corresponding  to  the  inter¬ 
face  roughness  (Figure  18)  seems  to  be  generating  the  worst  fading  environment.  This 
was  expected  since,  as  already  mentioned,  in  shallow-water  channels  the  interface  rough¬ 
ness  is  typically  the  dominant  scattering  factor. 

Figure  23  represents  the  worst-case  scenario  where  all  perturbations  are  combined 
and  Doppler  =  15  knots  «  8m/s)  is  added  to  account  for  the  spectral  broadening  in¬ 
duced  by  the  channel.  The  effect  of  Doppler  is  obvious  by  comparing  the  results  illus¬ 
trated  in  Figure  22  (all  perturbations  and  no  Doppler)  and  Figure  23.  It  is  obvious  that  the 
underwater  channel  is  far  more  challenging  for  communications  when  Doppler  is  in¬ 
cluded. 
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Parameter 

Value 

Remarks 

Main  Parameters 

Number  of  depth  points  saved 

3277 

Minimum  depth  saved 

Om 

Maximum  depth  saved 

200  m 

Number  of  range  steps  saved 

10 

Minimum  range  saved 

0  km 

Maximum  range  saved 

5  km 

Computational  range  step  size 

0.4  m 

Maximum  computational  depth 

500  m 

Number  of  computational  depth 
points 

16,384 

Radix-2  integer  required  for  FFT 

Reference  sound  speed 

1500  m/s 

Source  Data 

Source  depth 

75  m 

Center  frequency 

11.5  kHz 

Frequency  bandwidth 

2048  Hz 

Number  of  Frequencies 

1024 

Radix-2  integer  required  for  FFT 

Sound  Speed  Profile  Data 

Water  column  sound  speed 
(unperturbed) 

1500  m/s 

Isospeed 

Bathymetry  of  the  Water-Bottom  Interface 

Bottom  depth  (unperturbed) 

150  m 

Flat 

Acoustic  Parameters  of  the  Medium  below  the  Water-Bottom  Interface 

Bottom  sound  speed 

1600  m/s 

Sound  speed  gradient 

1/s 

Relative  density  (bottom/water) 

1.6 

Compressional  attenuation 

0.1  dB/km/Hz 

Shear  speed 

Om/s 

Shear  attenuation 

0  dB/m/kHz 

Acoustic  Properties  of  the  Deep  Layer 

Sub-bottom  sound  speed 

2000  m/s 

Sound  speed  gradient 

1/s 

Relative  density  (sub-bottom/water) 

3  glCYCl 

Compressional  attenuation 

0.25  dB/km/Hz 

Shear  speed 

Om/s 

Shear  attenuation 

0  dB/m/kHz 
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Parameter 

Value 

Remarks 

Perturbation  Parameters 

Source  speed 

15  knots  (8  m/s) 

Doppler  perturbation 

Bottom  interface  roughness  (rms) 

2  m 

Bottom  interface  perturbation 

Bottom  volume  sound  speed  vari¬ 
ability  (rms) 

15  m/s 

Bottom  volume  perturbation 
(with  corresponding  density  per¬ 
turbation,  as  described  in  text) 

Internal  wave  (linear)  sound 
speed  variability  (length  scales  & 
max  amplitudes) 

1500  m  &  1.5  m/s 

1200  m  &  1.5  m/s 

900  m  &  1.5  m/s 

600  m  &  1.5  m/s 

300  m  &  1.5  m/s 

First  part  of  IW  perturbation  (al¬ 
ways  combined  with  second  part 
oflW) 

Internal  wave  (soliton)  sound 
speed  variability  (ranges  &  max 
amplitudes) 

4800  m  &  12.5  m/s 

3800  m  &  9.3  m/s 

3000  m  &  6.9  m/s 

2400  m  &  5.1  m/s 

2000  m  &  3.8  m/s 

1800  m  &  2.8  m/s 

Second  part  of  IW  perturbation 
(always  combined  with  first  part 
oflW) 

Water  column  turbulence  (length 
scale  and  rms)  sound  speed  vari¬ 
ability 

100  m  &  2  m/s 

Turbulence  perturbation 

Table  1.  MMPE  Parameters. 
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Transmission  Loss  (dB  re  1m)  -  No  Perturbation 


1.05  1.1  1.15  1.2  1.25 

Freq  (Hz)  ^  ..a 


Normalized  Impulse  Response  -  No  Perturbation 


Travel  Time  (sec) 


Figure  17.  No  Perturbation;  TL  as  a  Funetion  of  Depth  and  Frequeney  (top  plot)  and 
Normalized  Impulse  Response  at  Depth  75  m  vs.  Travel  Time  (bottom  plot). 
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Transmission  Loss  (dB  re  1m)  •  Interface  Roughness  Perturbation 


Freq  (Hz)  ^ 


Normalized  Impulse  Response  -  Interface  Roughness  Perturbation 


Ill  'f|  II  N  J  'I  '  i 


3.1  3.2  3.3  3.4  3.5 

Travel  Time  (sec) 


3.6  3.7 


Figure  18.  Interface  Roughness  Perturbation;  TL  as  a  Function  of  Depth  and  Fre¬ 
quency  (top  plot)  and  Normalized  Impulse  Response  at  Depth  75  m  vs.  Travel  Time  (bot¬ 
tom  plot). 
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Transmission  Loss  (dB  re  1m)  •  Bottom  Volume  Perturbation 


Freq  (Hz)  ^ 


Normalized  Impulse  Response  -  Bottom  Volume  Perturbation 


3.1  3.2  3.3  3.4  3.5 

Travel  Time  (sec) 


3.6  3.7 


Figure  19.  Bottom  Volume  Perturbation;  TL  as  a  Function  of  Depth  and  Frequency 
(top  plot)  and  Normalized  Impulse  Response  at  Depth  75  m  vs.  Travel  Time  (bottom 

plot). 
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Transmission  Loss  (dB  re  1m)  •  Internal  Wave  Perturbation 


Travel  Time  (sec) 

Figure  20.  Internal  Wave  Perturbation:  TL  as  a  Function  of  Depth  and  Frequency 
(top  plot)  and  Normalized  Impulse  Response  at  Depth  75  m  vs.  Travel  Time  (bottom 

plot). 
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Normalized  Impulse  Response  -  Turbulence  Perturbation 
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Travel  Time  (sec) 


Figure  21.  Turbulence  Perturbation;  TL  as  a  Function  of  Depth  and  Frequency  (top 
plot)  and  Normalized  Impulse  Response  at  Depth  75  m  vs.  Travel  Time  (bottom  plot). 
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Normalized  Impulse  Response  -  All  Perturbations  (without  Doppler) 
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Figure  22.  All  Perturbations  Combined  (no  Doppler):  TL  as  a  Funetion  of  Depth  and 
Frequency  (top  plot)  and  Normalized  Impulse  Response  at  Depth  75  m  vs.  Travel  Time 

(bottom  plot). 
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Transmission  Loss  (dB  re  1m)  'All  Perturbations  (with  Doppler) 
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Figure  23.  All  Perturbations  Combined  (with  Doppler):  TL  as  a  Function  of  Depth 
and  Frequency  (top  plot),  and  Normalized  Impulse  Response  at  Depth  75  m  vs  Travel 

Time  (bottom  plot). 


Normalized  Impulse  Response  -  All  Perturbations  (with  Doppler) 
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_l _ 1 _ I_ 


This  chapter  presented  a  brief  description  of  the  MMPE  model  and  an  outline  of 
the  theory  behind  the  various  perturbations  modeled.  The  MMPE  model  was  used  to  pre¬ 
dict  the  channel  impulse  response  for  each  individual  perturbation  so  that  the  implemen¬ 
tation  of  the  DS-SS  scheme  with  the  CRV  updating  algorithm  could  be  evaluated  in 
Chapter  V.  In  the  following  chapter  the  CRV  decomposition  and  its  respective  updating 
algorithm  is  presented,  and  the  reason  that  the  CRV  can  be  used  as  an  alternative  to  the 
eigenvalue  decomposition  (EVD)  at  less  computational  cost  is  shown. 
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IV.  THE  CRV  DECOMPOSITION 


The  purpose  of  this  ehapter  is  to  deseribe  the  “Cross-produet  RV  deeomposition” 
[9].  The  eross-produet  RV  or  CRV  is  a  relatively  new  subspaee  traeking  algorithm  based 
on  the  URV  deeomposition  [29].  In  various  signal  proeessing  applieations,  but  mostly  in 
noise  suppression  teehniques,  it  is  often  needed  to  separate  the  signal  subspaee  from  the 
noise  subspaee.  The  most  powerful  tool  to  perform  this  task,  as  far  as  eomputational  ae- 
euraey  is  eoneerned,  is  the  singular  value  deeomposition  (SVD).  Although  the  SVD  is  a 
very  reliable  means  of  determining  the  numerieal  rank  of  a  eorrelation  matrix,  it  has  an 
important  drawbaek.  SVD  is  eomputationally  expensive,  espeeially  for  large  matriees, 
and  it  is  nonreeursive.  This  signifieant  disadvantage  ean  be  effeetively  addressed  by  the 
CRV. 

A.  CRV  DECOMPOSITION 

We  start  by  eonsidering  an  NxM  (N  >  M)  data  matrix  X  of  numerieal  rank  r  . 
Here,  numerieal  rank  means  the  number  of  singular  values  that  are  larger  than  a  eertain 
threshold.  Therefore  if  cTpCTj, . ,cr^  are  the  singular  values  of  X,  then 


CTj  >  CTj  > . threshold  »  ~ . «  cr^  «  0.  (4.1) 

By  definition  the  singular  values  cr,.  are  direetly  related  to  the  eigenvalues  A.  of  the  data 
eorrelation  matrix  X^X  sinee  A.  =  af .  Thus  the  number  of  eigenvalues  above  the  speei- 
fied  threshold  is  the  numerieal  rank  of  the  matrix.  In  signal  proeessing  if  the  gap  at  the 
threshold  is  well  determined,  we  say  that  the  r  largest  eigenvalues  are  assoeiated  to  the 
signal  subspaee  and  the  M  -r  smaller  eigenvalues  are  assoeiated  to  the  noise  subspaee. 
In  this  way,  eaeh  subspaee  is  spanned  by  the  eorresponding  eigenveetors.  The  CRV  de¬ 
eomposition  as  already  mentioned  is  based  upon  the  URV  algorithm  [29],  whieh  is  a 
rank-revealing  deeomposition  proposed  as  an  alternative  to  the  SVD.  The  URV  deeom¬ 
position  ean  provide  a  reliable  estimate  of  the  desired  subspaees  at  a  redueed  eomputa¬ 
tional  eost  and  may  be  written  as 


X  =  URV^ 


F 

G 


(4.2) 
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where  the  matriees  U,  V  are  unitary,  the  matriees  R,  E,  G  are  upper  triangular,  and  F  is 
a  matrix  of  very  small  numbers.  The  smallest  singular  value  (or  eigenvalue)  of  the£’  ma¬ 
trix  is  <j^  and  the  signal  subspaee  is  eoneentrated  in  this  matrix.  The  noise  subspaee  is 
eoneentrated  mainly  in  the  G  matrix. 

The  CRV  deeomposition  is  formulated  by  taking  the  produet  of  the  URV  and  its 
transpose.  If  as  before  X  is  an  NxM  data  matrix,  then  [9] 

0  =  X^X  =  [URV'^y  (URV")  =  VR^  [u^u)rV^.  (4.3) 

Considering  that  f/  is  a  unitary  matrix  {UU^  =  ij  and  by  setting 


W  =  R"R, 


(4.4) 


Equation  (4.3)  ean  be  written  as 


0  =  VWV^.  (4.5) 

Equation  (4.5)  is  a  eross-produet  RV  deeomposition.  By  eombining  (4.2)  and  (4.4)  the  W 
matrix  is  written  as 


~E 

F~ 

H 

~E 

F~ 

0  " 

~E 

F~ 

E^E 

E^F 

0 

G 

0 

G 

F" 

0 

G 

F^E 

F^F  +  G^G 

Einally  if  A  =  E^E,  B  =  E^ F  and  C  =  F^ F  +  G^ G  the  rank-revealing  CRV  deeompo¬ 


sition  is  formulated 


0  =  VWV^ =  V 


A  B 
B"  C 


(4.7) 


Sinee  the  smallest  singular  value  of  E  is  cr^,  the  smallest  eigenvalue  of  A  =  E^E  is 

/l^  =  cr^ .  Therefore  the  signal  subspaee  power  will  be  eoneentrated  mainly  in  matrix  A 

and  the  noise  subspaee  in  matrix  C. 

B,  CRV  UPDATING  ALGORITHM 

In  many  signal  proeessing  applieations,  the  rows  of  the  data  matrix  X  eorrespond 
to  samples  of  a  reeeived  signal  eorrupted  with  noise  and  the  eolumns  to  different  signals 
in  time.  If  the  signal  and  noise  subspaees  have  already  been  estimated  and  a  new  signal  is 


46 


received,  then  it  would  be  more  convenient  to  use  the  previously  estimated  noise  power 
to  approximate  the  new  subspaces  than  to  compute  everything  from  the  beginning.  This 
procedure  is  called  updating  [29]  and  can  be  implemented  recursively.  The  CRV  updat¬ 
ing  algorithm  described  next  is  a  simplified  version  of  the  algorithm  proposed  in  [9]. 

We  start  our  description  by  assuming  that  the  decomposition  0  =  VWV^  is  known 
and  the  updated  one  0  needs  to  be  computed.  The  updated  decomposition  is  given  by 

0  =  VWV"  =a0  +  xx" ,  (4.8) 

where  a  is  the  fading  (or  forgetting)  factor  with  values  between  0  and  1  and  x  the  newly 
received  M  x  1  data  vector.  In  our  implementation,  the  fading  factor  is  taken  to  be  1 .  Thus 
it  can  be  omitted  for  the  rest  of  the  algorithm  description.  The  next  step  is  to  define  the 
intermediate  vector  z  given  by 

z  =  F"x,  (4.9) 

so  that  equation  (4.8)  may  be  written  as, 

0  =  VWV'^  +Vz{Vzf  =VWV"  +  Vzz"V’'  =V  {w  +  zz'^Y"  =  V'FV‘' ,  (4.10) 

and  the  problem  is  transformed  into  updating  the  matrix  ^  =  W  +  zz"  .  If  the  numerical 
rank  of  W  according  to  a  user  specified  threshold  is  r  then  IP  is  of  the  form  (example  for 
r  =  3): 

X  X  X  e  e 

X  X  X  e  e 

X  X  X  e  e 

T  =  e  e  e  e  e 

e  e  e  e  e 


\e  e  e  e  e  ■■■  e 
w 

where  the  upper  left  rxr  elements  of  W,  denoted  with  x ’s,  will  have  very  large  values 
compared  to  the  rest  (^M  -r^  of  the  elements,  denoted  with  e’s.  Now  that  the 

matrix  IP  is  defined,  the  next  step  in  the  algorithm  is  to  find  a  Householder  matrix  P  and 
apply  a  transformation  to  keep  the  first  r  + 1  elements  of  z  by  clearing  the  last  M  -r-\ 
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elements.  Householder  matriees  P  (or  elementary  refleetors)  [31]  have  the  important  abil¬ 
ity  to  zero  speeified  entries  in  a  matrix  or  veetor  and  are  defined  as 


with  /  the  identity  matrix, 


P  =  I-2 


vv 

T  • 

V  V 


and 


so  that 


V  =  z± 


1 

0 


Pz  =  ±  z 


(4.12) 


(4.13) 


(4.14) 


(4.15) 


then 

HL 

=  5,  v  = 

[7,  1,  4,  2]^  and 

~2 

1 

4 

2 

1 

1 

-34/7 

4/7 

2/7 

5 

4 

4/7 

-19/7 

8/7 

2 

2/7 

8/7 

-31/7 

sueh  that 

Pz  =  -5ei=[-5,  0,  0,  of. 

After  eomputing  the  (M-r)x(M-r)  Householder  matrix  P,  needed  to  elear  the  last 
M  -r-\  entries  of  veetor  z,  we  form  the  M  xM  bloek  diagonal  matrix  P  as 

L  0 


P  = 


x(M-r) 

0  P 

^(M-r)xr  (M-r)x[M-r) 


(4.16) 


With  the  help  of  equation  (4.16)  the  updated  matrix  P  is  given  by 

P  =PPP^  =p(w  +  zz^)p^ , 


(4.17) 
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and  is  of  the  form 


<//  = 


X 

X 

X 

X 

e  ••• 

e 

X 

X 

X 

X 

e  ••• 

e 

X 

X 

X 

X 

e  ••• 

e 

X 

X 

X 

X 

e  ••• 

e 

e 

e 

e 

e 

e  ••• 

e 

e 

e 

e 

e 

e  ••• 

e 

(4.18) 


The  upper  left  (r +  l)x(r +  l)  submatrix  of  *F  denoted  with  x  5  is  unknown  and,  in  order 


to  proeeed  with  the  next  eyele  of  the  iterative  algorithm,  the  rank  f  of  'F  needs  to  be  es¬ 
timated.  This  is  done  by  eomputing  the  eigenvalues  of  the  (r-l-l)x(r +  l)  bloek  of  x’s 

and  eomparing  the  two  smallest  ones  by  putting  a  simple  deeision  threshold.  Speoifieally, 
if  /l^^j  >0.1/1^  then  we  determine  that  the  rank  has  ehanged  and  r  =  r  +  \.  On  the  other 


hand  if  /l^^j  <0.1/1^  the  rank  of  the  updated  matrix  "F  remains  unehanged  and  r  =  r. 

The  way  that  the  rank  changes  are  determined  in  this  implementation  of  the  CRV  updat¬ 
ing  algorithm  is  not  the  most  efficient  one  as  far  as  computations  are  concerned.  For  “r  ” 
small  (say  2  or  3)  there  are  closed  formulas  that  can  be  applied  for  the  computation  of  the 
eigenvalues,  as  the  Lanczos  bidiagonilization  algorithm  [10].  These  formulas  offer  the 
same  accuracy  at  less  computational  cost  but  are  not  examined  in  this  research.  All  the 
steps  described  above  are  repeated  each  time  a  new  signal  (vector  x)  is  received.  At  the 
end  of  each  cycle,  the  correlation  matrix  0  is  computed  from  Equation  (4.10)  and  is 
given  by 


0  =  V4'V^ , 


(4.19) 


where 


V  =  VP.  (4.20) 

A  graphical  representation  of  the  CRV  updating  algorithm  is  illustrated  in  the  fol¬ 
lowing  flowchart  (Figure  24). 
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Figure  24.  CRV  Updating  Algorithm. 
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C.  NUMERICAL  EXAMPLE 

In  this  section  a  simple  numerical  example  (done  with  MATLAB)  is  presented  to 
demonstrate  the  accuracy  of  the  correlation  matrix  estimate  produced  with  the  CRV  up¬ 
dating  algorithm  illustrated  in  Figure  24.  We  start  by  creating  a  10x10  data  matrix  X 
Each  column  of  X  represents  a  random  signal  x. ,  i  =  1,2,  ....,10  of  ones  and  minus  ones 
corrupted  with  additive  white  noise: 


-0.315 

-0.781 

-2.837 

1.506 

-2.731 

-0.006 

-0.223 

0.654 

-0.199 

-1.379 

0.698 

-1.66 

-0.502 

-2.028 

1.478 

0.306 

-0.696 

1.986 

1.88 

-1.242 

-1.599 

2.351 

-0.18 

-0.9 

0.552 

2.199 

-1.73 

1.643 

0.775 

-0.165 

-0.219 

1.136 

0.765 

0.833 

1.386 

-3.577 

-2.143 

1.919 

-0.703 

1.756 

1.602 

-0.833 

-0.631 

1.689 

1.053 

-1.086 

-0.413 

-2.248 

-1.521 

2.164 

-0.057 

-0.705 

-1.317 

0.883 

0.513 

-0.613 

0.408 

-0.842 

0.844 

-2.023 

-0.023 

1.276 

0.203 

-0.674 

1.244 

0.138 

1.518 

1.788 

0.901 

0.701 

0.932 

1.394 

-0.309 

0.904 

1.718 

-2.23 

-2.492 

0.422 

-0.002 

-1.494 

1.081 

-1.098 

0.957 

1.031 

-0.846 

1.641 

0.913 

-0.472 

-0.565 

-0.827 

-2.767 

1.176 

1.324 

-1.613 

-0.866 

0.095 

0.987 

2.671 

0.974 

1.354 

Xj 

Xj 

Xj 

XiQ 

The  correlation  matrix  estimate  is  computed  for  the  signals  Xj  and  Xj  by  going 

through  all  the  steps  of  the  first  two  cycles  of  the  algorithm. 

1,  First  Cycle 

The  algorithm  is  initialized  by  setting  W  =  0  andV  =  I .  Both  matrices  are  of  di¬ 
mension  10x10  and  the  rank  r  of  IF  is  zero.  Next  we  need  to  compute  the  z  vector  given 
by  Equation  (4.9): 


z  =  F^X[  =  /jqXj 


-0.315 

0.698 

-1.599 

-0.219 

1.602 

-0.057 

-0.023 

0.932 

1.081 

-2.767 


(4.22) 
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where  x^  is  the  first  eolumn  of  the  data  matrix  X  From  Equation  (4.16)  and  the  fact  that 

r  =  0  we  get  P  =  P ,  where  P  is  the  Householder  matrix  needed  to  clear  (set  to  zero)  all 
but  the  first  element  of  z.  The  Householder  matrix  P  is  computed  with  the  help  of  the  ex¬ 
pressions  (4.12)  through  (4.14).  Now  that  P  is  defined  the  updated  matrix  P  given  by 
(4.17)  is 

15.463  0  •••  0 
0  0  •••  0 

0  0  •••  0 


.  (4.23) 

10x10 


.  X  w=0 

P  =ppp^  =ppp^  =p{w  +  zz^)p"  =Pzz"P^  = 


and  the  updated  matrix  V  by  (4.20): 


-0.08 

0.177 

-0.406 

-0.055 

0.407 

-0.014 

-0.006 

0.237 

0.275 

-0.703 

0.177 

0.97 

0.066 

0.009 

-0.067 

0.002 

0.001 

-0.038 

-0.045 

0.115 

-0.406 

0.066 

0.846 

-0.021 

0.153 

-0.005 

-0.002 

0.089 

0.103 

-0.264 

-0.055 

0.009 

-0.021 

0.997 

0.021 

-0.0007 

-0.0003 

0.012 

0.014 

-0.036 

0.407 

-0.067 

0.153 

0.021 

0.846 

0.005 

0.002 

-0.089 

-0.103 

0.265 

-0.014 

0.002 

-0.005 

-0.0007 

0.005 

0.999 

0.000 

0.003 

0.003 

-0.009 

-0.006 

0.001 

-0.002 

-0.0003 

0.002 

0.000 

0.999 

0.001 

0.001 

-0.003 

0.237 

-0.038 

0.089 

0.012 

-0.089 

0.003 

0.001 

0.947 

-0.06 

0.154 

0.275 

-0.045 

0.103 

0.014 

-0.103 

0.003 

0.001 

-0.06 

0.929 

0.179 

-0.703 

0.115 

-0.264 

-0.036 

0.265 

-0.009 

-0.003 

0.154 

0.179 

0.541 

(4.24) 


The  last  step  is  to  compute  the  correlation  matrix  estimate  O  (4.19): 


0.099 

-0.22 

0.505 

0.069 

-0.506 

-0.22 

0.487 

-1.117 

-0.153 

1.119 

0.505 

-1.117 

2.558 

0.351 

-2.564 

0.069 

-0.153 

0.351 

0.048 

-0.352 

-0.506 

1.119 

-2.564 

-0.352 

2.569 

0.018 

-0.039 

0.091 

0.012 

-0.091 

0.007 

-0.016 

0.038 

0.005 

-0.038 

-0.294 

0.651 

-1.491 

-0.205 

1.494 

-0.341 

0.755 

-1.73 

-0.237 

1.734 

0.873 

-1.932 

4.426 

0.608 

^.435 

0.018 

0.007 

-0.295 

-0.341 

0.873 

-0.039 

-0.016 

0.651 

0.755 

-1.932 

0.091 

0.038 

-1.491 

-1.73 

4.426 

0.012 

0.005 

-0.205 

-0.237 

0.608 

-0.091 

-0.038 

1.494 

1.734 

^.435 

0.003 

0.001 

-0.053 

-0.061 

0.158 

0.001 

0.0005 

-0.022 

-0.025 

0.066 

-0.053 

-0.022 

0.868 

1.008 

-2.579 

-0.061 

-0.025 

1.008 

1.17 

-2.993 

0.158 

0.066 

-2.579 

-2.993 

7.656 

(4.25) 
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The  rank  r  of  'F  is  determined  (only  in  the  first  cycle  of  the  algorithm)  with  the  Matlab 
function  “rank,”  which  provides  an  estimate  of  the  number  of  linearly  independent  rows 
or  columns  of  a  matrix.  In  this  example  the  rank  r  of  IP  is  one. 

2,  Second  Cycle 

The  CRV  updating  algorithm  is  a  recursive  one;  thus  the  updated  matrices  V  and 
W  =  W  computed  in  the  first  cycle  are  used  to  initialize  the  second  cycle  by  setting 
W  =  W ,V  =  V  and  r  =  \.  The  z  vector  is  computed  from  (4.9): 

-2.388 
-1.396 

1.746 
1.053 
-0.226 
-0.726 
1.267 

1.747 
-0.689 
0.129 

where  Xj  is  the  second  column  of  X.  The  Householder  matrix  P  needed  to  zero  all  except 
the  first  two(r +  1)  elements  of  z  is  determined  from  Equations  (4.12)  through  (4.14).  For 
computational  efficiency  P  is  determined  by  operating  on  the  last  9  (10 -r)  elements  of  z 
and  is  of  the  dimension  9x9: 

■-0.405  0.507  0.306  -0.065  -0.211  0.368  0.507  -0.200  0.037 

0.507  0.816  -0.110  0.023  0.076  -0.132  -0.183  0.072  -0.013 

0.306  -0.110  0.933  0.014  0.046  -0.080  -0.110  0.043  -0.008 

-0.065  0.023  0.014  0.996  -0.009  0.017  0.023  -0.009  0.001 

P=  -0.211  0.076  0.046  -0.009  0.968  0.055  0.076  -0.030  0.005 

0.368  -0.132  -0.080  0.017  0.055  0.903  -0.133  0.052  -0.009 

0.507  -0.183  -0.110  0.023  0.076  -0.133  0.816  0.072  -0.013 

-0.200  0.072  0.043  -0.009  -0.030  0.052  0.072  0.971  0.005 

0.037  -0.013  -0.008  0.001  0.005  -0.009  -0.013  0.005  0.998 
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The  block  diagonal  matrix  P  is  determined  from  Equation  (4.16)  as 


Ir 

O 

1 

“A 

*^1x9 

^(M-r)xr 

1 

1 

1 _ 

1 

o 

so 

X 

^9x9  _ 

where  P  is  the  Householder  matrix  given  by  (4.27).  The  resulting  matrix  is 


(4.28) 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.405 

0.507 

0.306 

-0.065 

-0.211 

0.368 

0.507 

-0.200 

0.037 

0 

0.507 

0.816 

-0.110 

0.023 

0.076 

-0.132 

-0.183 

0.072 

-0.013 

0 

0.306 

-0.110 

0.933 

0.014 

0.046 

-0.080 

-0.110 

0.043 

-0.008 

0 

-0.065 

0.023 

0.014 

0.996 

-0.009 

0.017 

0.023 

-0.009 

0.001 

0 

-0.211 

0.076 

0.046 

-0.009 

0.968 

0.055 

0.076 

-0.030 

0.005 

0 

0.368 

-0.132 

-0.080 

0.017 

0.055 

0.903 

-0.133 

0.052 

-0.009 

0 

0.507 

-0.183 

-0.110 

0.023 

0.076 

-0.133 

0.816 

0.072 

-0.013 

0 

-0.200 

0.072 

0.043 

-0.009 

-0.030 

0.052 

0.072 

0.971 

0.005 

0 

0.037 

-0.013 

-0.008 

0.001 

0.005 

-0.009 

-0.013 

0.005 

0.998 

and  the  updated  matrix  W  =  PWP^  =  P[w  +  zz^^P^  produces  the  following: 


“  21.17 

-8.219 

0 

...  o' 

-8.219 

11.838 

0 

...  0 

0 

0 

0 

...  0 

9 

(4.30) 

0 

0 

0 

...  0 

10x10 

where  the  entries  of  W  ,  denoted  with  zeros,  are  very  small  numbers (^<  10  .  The  eigen¬ 

values  of  the  upper  left  2x2(r  +  lxr  +  l)  block  of  W  are  \  =25.955  and  =7.052 


and,  since /Ij  >  0.1/l[,  we  determine  that  the  rank  has  changed  by  one  and  is 

now r  =  r-l-l  =  2.  It  must  be  noted  here  that  the  decision  threshold  depends  on  the  appli¬ 
cation,  and  there  is  no  specific  method  known  so  far  that  could  lead  to  an  optimum  choice 
[9].  The  signal  subspace  power  is  concentrated  in  the  upper  left  2x2  block  of  W  whereas 
the  noise  subspace  is  in  the  lower  right  8x8  block.  Finally,  the  correlation  matrix  esti¬ 
mate  0  calculated  from  Equations  (4.19)  and  (4.20)  is  the  following: 
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0.710 

1.076 

-1.331 

-0.818 

0.144 

0.569 

-0.989 

-1.383 

0.516 

-0.045 

1.076 

3.244 

-5.021 

-2.040 

2.503 

1.131 

-2.136 

-1.664 

2.579 

-3.886 

-1.331 

-5.021 

8.087 

3.023 

^.523 

-1.566 

3.039 

1.788 

^.313 

7.192 

-0.818 

-2040 

3.023 

1.339 

-1.299 

-0.788 

1.455 

1.379 

-1.486 

1.945 

0.144 

2.503 

^.523 

-1.299 

3.263 

0.495 

-1.102 

0.332 

2.649 

-5.415 

0.569 

1.131 

-1.566 

-0.788 

0.495 

0.500 

-0.898 

-1.036 

0.712 

-0.671 

-0.989 

-2.136 

3.039 

1.455 

-1.102 

-0.898 

1.630 

1.757 

-1.428 

1.567 

-1.383 

-1.664 

1.788 

1.379 

0.332 

-1.036 

1.757 

2.813 

-0.523 

-0.938 

0.516 

2.579 

^.313 

-1.486 

2.649 

0.712 

-1.428 

-0.523 

2.377 

^.285 

-0.045 

-3.886 

7.192 

1.945 

-5.415 

-0.671 

1.567 

-0.938 

^.285 

9.040 

The  exact  correlation  matrices  for  the  first  two  signals  x^  and  Xj  are  determined 
with  the  following  expressions: 

(4.32) 

and 

•  (4.33) 

After  carrying  out  the  calculations,  the  results  from  Equations  (4.32)  and  (4.33) 
are  exactly  the  same  with  the  correlation  matrix  estimates  (4.25)  and  (4.31)  computed 
with  the  CRV  updating  algorithm.  The  eigenvalues  of  the  exact  correlation  matrix  and  the 
eigenvalues  of  the  estimated  one  for  the  l'^*,  3’^‘*,  6**'  and  10*  iteration  are  illustrated  in 
Figure  25.  By  looking  at  the  figure,  we  conclude  that  the  estimate  produced  by  the  CRV 
algorithm  is  very  accurate. 
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Eigenvalues  of  Correlation  Matrix  Estimate  Computed  With  the  CRV  Algorithm 
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Eigenvalues  of  the  Exact  Correlation  Matrix 
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Figure  25.  Values  of  the  n-th  Eigenvalue  of  the  Estimated  and  the  Exact  Correlation 

Matrix. 

In  this  chapter,  the  CRV  updating  algorithm  [9]  for  subspace  tracking  was  de¬ 
scribed  and  a  simple  numerical  example  was  presented  so  that  its  accuracy  could  be 
evaluated.  In  the  next  chapter  of  this  thesis  the  performance  of  the  algorithm  is  examined 
by  implementing  it  in  a  previously  proposed  communication  scheme  [2]  under  varying 
channel  conditions. 
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V.  IMPLEMENTATION  AND  PERFORMANCE  OF  THE  CRV 

UPDATING  ALGORITHM 


This  chapter  describes  the  implementation  of  the  CRV  updating  algorithm  in  a 
previously  proposed  direct  sequence  spread-spectrum  (DS-SS)  communication  scheme 
[2].  The  performance  of  the  algorithm  is  evaluated  under  different  channel  conditions  and 
design  parameters  (packet  size,  bit  rate).  The  theory  relevant  to  spread-spectrum  (SS) 
signaling  is  presented  first  and  the  receiver  structure  is  described  next.  The  performance 
results  are  presented  last  in  terms  of  Bit-Error-Rate  (BER)  vs.  Bit  Energy  per  Noise 
Power  as  it  is  usually  done  in  all  communication  systems. 

A.  DIRECT  SEQUENCE  SPREAD  SPECTRUM  (DS-SS) 

Spread-spectrum  (SS)  techniques  are  preferred  in  military  communications  be¬ 
cause  of  their  low  probability  of  intercept  (EPI),  low  probability  of  detection  (EPD)  and 
their  anti-jam  capability  [8].  In  SS  applications  the  resulting  bandwidth  of  the  transmitted 
signal  is  much  greater  than  the  information  bandwidth  and  the  transmitted  signal  appears 
to  an  unauthorized  receiver  as  noise.  The  two  most  commonly  used  SS  signaling  schemes 
are  frequency  hopping  (PH)  and  direct  sequence  (DS).  Only  DS  is  described  in  this  the¬ 
sis.  The  block  diagram  of  a  simple  DS-SS  system  is  illustrated  in  Figure  26. 


Figure  26.  A  Simple  DS-SS  System. 


The  binary  message  signal  is  spread  at  baseband  by  a  locally  generated  pseudo¬ 
noise  code  (PN-code).  The  PN  sequence  c{t) ,  also  known  as  the  chipping  sequence,  is  a 
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binary  signal  (c(0  =  ±l)  that  appears  random  and  is  produced  at  much  higher  frequency 
than  the  data  that  is  to  be  transmitted.  That  is, 

K»R„  (4.34) 


where  is  the  chip  rate  in  chips-per-second  [cps]  and  is  the  information  data  rate  in 
bits-per-second  [bps].  Since  the  chip  duration  7)  is  much  less  than  the  bit  dura¬ 

tion  .  The  processing  gain  (or  multiplicity  factor)  N  (always  an  integer)  in  chips-per-bit 
[cpb]  for  a  DS  system  is 


(4.35) 


Therefore  the  message  spectrum  is  spread  N  times  before  the  message  is  modulated  and 
transmitted.  The  wide  bandwidth  provided  by  the  PN  code  allows  the  signal  to  be  hidden 
in  noise  since  its  power  is  dropped  below  the  noise  threshold.  In  this  way  SS  signals  are 
very  difficult  to  detect  (LPD).  The  process  of  spreading  by  multiplication  with  a  PN  se¬ 


quence  is  shown  in  Figure  27  for  the  case  of  A^=  2  (7)  =  0.57],)  . 


Figure  27.  Example  of  Spreading  in  a  DS-SS  System. 
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There  are  many  different  PN  codes  used  in  DS-SS  communications  like  the 
maximal  length  codes,  the  Hadamard- Walsh  codes,  the  Gold  codes,  the  Kasami  codes, 
and  others.  In  the  implementation  of  this  thesis  research,  the  spreading  of  the  data  signal 
is  done  by  Gold  codes,  which  are  produced  with  the  Gold  code  generators  of 
SIMULINK.  Gold  codes  have  very  good  cross-correlation  properties  (the  cross¬ 
correlation  of  two  different  Gold  codes  is  very  small)  and  are  preferred  when  multiple 
users  must  share  the  same  frequency  spectrum  (each  user  is  assigned  a  different  code).  As 
shown  in  Figure  21,  the  exact  replica  of  the  PN  code  used  at  the  transmitter  is  needed  in 
order  for  the  receiver  to  decode  the  transmitted  signal.  The  decoding  is  done  by  a  simple 
multiplication  (de-spreading)  of  the  received  signal  with  the  chipping  code 
sincec^(t)  =  l  (c(t)  =  ±l).  The  fact  that  the  receiver  must  be  able  to  generate  the  same 

spreading  sequences  in  a  deterministic  way,  just  as  the  transmitter  does,  provides  a  cer¬ 
tain  degree  of  communication  security.  The  transmitted  signal  if  detected  by  a  non¬ 
intended  receiver  cannot  be  intercepted  (decoded)  without  prior  knowledge  of  the  spe¬ 
cific  PN  sequence  used  for  the  spreading  of  the  information  signal  (LPI). 

The  third  important  advantage  in  using  SS  signaling  is  the  ability  to  resist  hostile 
interference  (jamming).  A  useful  parameter  in  specifying  the  performance  of  a  spread- 
spectrum  system  in  the  presence  of  jamming  is  the  processing  gain  N.  Since  the  message 
signal  is  spread  before  transmission,  the  interference  bandwidth  occupies  only  a  small 
portion  of  the  SS  signal  bandwidth,  as  shown  in  Figure  28  (left  plot).  At  the  receiver,  af¬ 
ter  de-spreading  the  received  signal  by  multiplying  with  the  same  PN  code  (the  one  used 
at  the  transmitter),  the  spectrum  of  Figure  28  (right  plot)  is  produced.  Most  of  the  origi¬ 
nal  interference  energy  is  eliminated  by  spreading  (except  for  the  energy  corresponding 
to  the  interference  spectrum  overlapping  with  the  signal  spectrum)  and  the  original  signal 
can  be  recovered.  It  is  also  clear  from  the  same  figure  that  the  greater  the  processing  gain 
N  of  the  SS  system,  the  better  is  the  performance  of  the  system  against  interference  (bar¬ 
rage  noise  jamming)  [18]. 
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Spectral  Density 


Spectral  Density 


Figure  28.  Spectrum  of  the  Transmitted  SS  Signal  (left  plot)  and  of  the  Received  Sig¬ 
nal  (right  plot)  after  De-spreading  [From  Ref.  18]. 


B,  PREVIOUSLY  PROPOSED  DS-SS  SCHEME 

The  channel  equalization  proposed  by  [2]  is  based  on  multi-rate  signal  processing 
and  subspace  decomposition.  The  analysis  is  done  in  the  discrete  time  domain  and  the 
transmitted  signal  is  assumed  to  be  a  binary  sequence  (+1,-1),  which  is  then  chipped  by  a 
short  PN  code  and  transmitted  through  the  channel.  Also,  error  correction  coding  and 
modulation  were  not  examined  and  are  not  implemented  in  this  thesis.  Thus  the 
worst-case  scenario  is  simulated  since  both  processes  (error  correction  coding  and  modu¬ 
lation)  can  only  improve  the  equalization  mechanism.  The  block  diagram  of  the  simpli¬ 
fied  DS-SS  system  is  shown  in  Figure  29. 


Channel  &  Interference 


Figure  29.  Simplified  DS-SS  System  [After  Ref.  2]. 


The  information  sequence  m[n]  is  first  up-sampled  by  N  (addA-1  zeros  be¬ 
tween  bits)  in  order  to  match  the  chipping  rate  and  then  is  spread  by  convolution  with 
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the  chipping  sequence  c\h\ .  Spreading  can  also  be  done  by  multiplying  the  up-sampled 
(where  each  bit  is  repeated  N  times)  data  sequence  m[n\  with  c{n\ .  The  spread  signal 
s\n\  is  then  convolved  with  the  channel  impulse  response  h\n\  and  noise  u{n]  is  added 
(white  or  colored)  to  account  for  the  random  interference.  Subsequently,  the  received  se¬ 
quence  is  de-spread  at  the  receiver  by  convolution  with  a  time-reversed  replica  c[-n]  of 
the  PN-code  c\n\  and  then  is  down-sampled  by  N.  In  order  to  effectively  remove  the 
“pseudo-randomization”  imposed  on  the  data  sequence  by  the  spreading  code,  a  convolu¬ 
tion  between  c{n\  and  c{-n\  must  yield  a  delta  function. 


c[n]*c[-n]  =  d\n\. 


(4.36) 


If  c{n\  =  [cq,  Cj, . ,  c^_i],  where  Cq,  Cj, . ,  c^_j  are  the  chips  of  the  PN-code,  then 


c[-n\  =  [: 


(4.37) 


For  example,  when  a  Gold  code  c\n\  of  length  125  is  used  and  c[-n]  is  given  by  Equa¬ 
tion  (4.37),  the  result  of  (4.36)  is  the  impulse  shown  in  Figure  30. 
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Figure  30.  The  Result  of  c[n]  *c[-n]  When  Gold  Codes  are  Used. 

The  equivalent  state-space  structure  of  the  communication  system  shown  in  Fig¬ 
ure  29  is  illustrated  in  Figure  3 1 . 
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A  A  A  A  A  AAA 

Transmitter  Channel  +  Interference  Receiver  Equalizer  ; 


Figure  31.  State-Space  Structure  of  the  DS-SS  System  [After  Ref  2], 

The  spreading  of  the  information  sequence  m[n]  is  accomplished  by  filtering  the 

up-sampled  (by  N)  data  sequence  by  a  A-th  order  linear  filter  defined  (in  the  z-domain)  as 

C(z)  =  c,+c,  z-'  +  Cj  Z-"  + . +  (4.38) 

The  spread  signal  s\n\  is  then  filtered  by  a  second  A-th  order  linear  filter  (channel  im¬ 
pulse  response), 

H{z)  =  Hq  +1\z^^  -t . + (4.39) 

to  account  for  the  corruption  caused  by  the  underwater  acoustic  channel.  At  the  receiver 
the  received  sequence  is  de-spread  by  multiplication  (and  not  convolution)  with  a  replica 
of  the  PN-code  used  at  the  transmitter.  In  practice  the  length  of  the  received  sequence 
does  not  change  after  de-spreading.  The  summation  step  implied  by  the  convolution  ap¬ 
proach  (c[n]*c[-n])  is  performed  later  at  the  equalizer  when  the  estimate  m\n\  of  m\n\ 

is  computed.  At  the  equalizer  (see  Figure  31),  the  de-spread  sequence  is  forward  delayed 
( z  elements  in  the  state-space  structure)  and  down-sampled  by  A  .  The  result  of  this  proc¬ 
ess  is  A  different  versions  of  the  original  information  sequence  (corrupted  by  the  fading 
channel)  denoted  asyo[n],  y\n\, . ,  .  For  a  better  understanding  of  how  the 

structure  of  Figure  31  works,  simple  examples  are  presented  in  Figure  32  (transmitter 

section)  and  Figure  33  (channel  and  receiver  section). 
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Figure  32.  Up-sampling  and  Spreading  by  Multiplication  with  the  PN-code. 

The  data  sequence  m[«]  consists  of  5  bits  [l,  1,  -1,  1,  -l]  and  the  processing  gain 
is  assumed  to  be  three (A^  =  3) .  Since =  3,  m{n\  is  up-sampled  by  3  (each  bit  is  re¬ 
peated  three  times)  and  then  bit-by-bit  multiplied  with  the  PN-code  c\n\  to  yield  the 
spread  sequence  s{n\ .  Subsequently,  white  noise  w[n]  is  added  to  account  for  the  chan¬ 
nel  corruption.  The  sequence  s\n\  +  w[n]  is  then  de-spread  at  the  receiver  by  multiplica¬ 
tion  with  a  replica  of  c\n\  (Figure  33,  third  plot  from  top).  Finally,  the  de-spread  se¬ 
quence  is  forward  delayed  and  down-sampled  by  N,  so  that  three  different  ver¬ 
sions  y\n\,  y2{n\  of  the  original  data  sequence  are  produced  at  the  equalizer  input. 

The  three  sequences  y^in],  y\n\,  y2{n\  are  shown  in  Figure  33  (first  plot  from  bottom)  in 
different  colors  (blue,  red,  and  green  respectively). 
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Channel  and  Receiver  Section 


Figure  33.  Channel  Corruption,  De-spreading  and  Different  Versions  of  the  Original 

Data  Sequence. 

The  vectors  y\n\, . ,  (see  Figure  31)  form  the  data  matrix 

(4.40) 

which  is  then  used  to  determine  the  signal  and  noise  subspace  (the  most  important  part  in 
the  equalization  process).  The  estimated  sequence  m\n\  is  given  by 

m{n-\^f-y{nl  (4.41) 

where  /  is  the  optimum  filter  that  maximizes  the  signal-to-noise  ratio  at  the  receiver. 
The  channel  equalization  is  performed  with  two  different  blind  equalization  algorithms, 
depending  on  whether  the  synchronization  between  the  receiver  and  the  transmitter  was 


y[n\  = 


-  To[«] 

-  y\n\ 

-  yN-M 
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achieved  or  not.  The  two  algorithms  are  based  on  matehed  filter  theory  and  are  briefly 
deseribed  in  the  following  subsections.  A  detailed  derivation  and  deseription  can  be 
found  in  Referenee  [2]. 

1.  Demodulation  with  Synchronization 

This  is  the  ease  where  the  replica  of  the  PN-code  used  at  the  transmitter  is  com¬ 
pletely  aligned  (in  time)  with  the  reeeived  sequenee.  Thus  the  de-spreading  is  performed 
correetly. 


a.  Multipath  with  Additive  White  Gaussian  Noise  (A  WGN) 

The  received  sequence  y[n]  is  equal  to  s[n]*  h[n]  + M{n]  where  s[n]  is 
the  spread  by  the  PN-eode  information  sequence,  w[n]  is  the  AWGN,  and  h[n]  is  the 
channel  impulse  response.  It  must  be  noted  here  that  the  maximum  length  of  the  impulse 
response  h[n] ,  according  to  the  theoretical  approach  described  in  the  previous  section,  is 
N,  where  N  is  the  processing  gain  (see  Figure  31).  If  m[n\  is  an  estimate  of  the  informa¬ 
tion  sequence  m[n]  then 

m{n-\  =  f-y{n\,  (4.42) 

where  /  is  the  eigenvector  corresponding  to  the  maximum  eigenvalue  of  the  data  eorre- 
lation  matrix  given  by 

Ryy=E[y[n-\/[n-\].  (4.43) 

b.  Multipath  with  Additive  Colored  Noise  (ACN) 

The  reeeived  sequenee  y{n\  is  equal  to  s{n\  *  h{n\  +  v{n\  where  v[n]  is  the 
additive  eolored  noise.  In  order  to  be  able  to  follow  the  approach  mentioned  above  in 
(4.42),  a  simultaneous  diagonalization  of  the  correlation  matrices  of  the  signal  and 

the  noise  R^  is  needed.  This  is  done  with  help  of  the  Mahalanobis  transformation,  de¬ 
fined  as  [32] 

y\n\  =  M-^'^y{nl  (4.44) 
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where  M  is  the  transformation  matrix  given  by 

M-^'^=EX^"eI=K:'\  (4.45) 

and  are  the  matrices  produced  from  the  eigenvalue  decomposition  of 

=£'|v[n]v^[n]| .  Specifically  E^  is  a  matrix  with  columns  of  the  eigenvectors  of 
and  is  a  diagonal  matrix  whose  elements  are  the  corresponding  eigenvalues.  The  cor¬ 
relation  matrixi?^of  y'\n\  (4.44)  is  computed  from 

R'yy  =E[y\n\{y\n-\)^]  =  (4.46) 

and,  since  =  R^l'^ ,  from 


R'„=R:1'^R^X'"-  (4-4V) 

Because  of  this  whitening  transformation  in  the  new  coordinate  system  the 
colored  noise  is  considered  white  and  the  estimate  m\n\  of  the  data  signal  m[n\  is  given 
by  (4.42) 

=  y\nl  (4.48) 

where  /'  is  the  eigenvector  corresponding  to  the  maximum  eigenvalue  of  the  autocorre¬ 
lation  matrix  R'^ . 

2.  Demodulation  without  Synchronization 

The  underwater  communication  channel  is  very  challenging  and  synchronization 
is  not  always  possible.  The  slightest  misalignment  in  the  PN-code  can  incorrectly  de¬ 
spread  the  received  sequence.  When  synchronization  is  achieved,  the  data  correlation  ma¬ 
trix  i?^(4.43)  has  only  one  dominant  eigenvalue  and  the  signal  subspace  is  easily  dis¬ 
criminated  from  the  noise  subspace.  However,  this  is  not  the  case  when  synchronization 
is  not  possible.  Suppose  that,  at  a  given  time,  y[n]  contains  two  symbols,  the  current 

symbol  m[n]  and  a  portion  of  the  previous  symbol  m[n-l\ ,  i.e. 


y[n]  =  gom[n]  +  g^m[n  - 1]  =  [go 


m[n] 

m[n-l\ 


(4.49) 
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where  g^,  gj  are  veetors  depending  on  the  autoeorrelation  of  the  PN-eode  and  the  eom- 
munieation  channel.  In  this  case  y[n]  spans  a  two-dimensional  subspace  defined  by  the 
vectors  gg  and  gj  and  there  are  two  dominant  eigenvalues  instead  of  one  in  the  data  cor¬ 
relation  matrix  R^.  If  A^[n],  A^[n]  are  the  two  dominant  eigenvalues  and  ^  their  re¬ 
spective  eigenvectors,  y[n]  is  given  by: 


y[n]  =  ^^n]  + 


(4.50) 


Since  the  eigenvectors  are  orthonormal  =  1  and^^^  =  0  ),  it  is  true  that 


and 


Ai[n]  =  ^^  -yin] 


A^[n]  =  §2^  ■  y[n]. 


By  combining  Equations  (4.51)  and  (4.52)  in  matrix  form,  we  get 


A,[n] 

_^[«r 

_eJ_ 

y[n\<^A[n\  = 


y[nl 


Substituting  (4.49)  into  (4.53)  yields 


A[n\  = 

If  we  solve  Equation  (4.54)  for 


•_£o  £i_' 


m[n\ 

m\n-\\ 


m[n] 

m\n-\\ 


.  we  get 


m\n\ 

m\n-\\ 


J  VL^ 


•go  gi_ 


•l[n]  <:^> 


) 

T  j  T 


m[n\ 

m[n-l\ 


ki 


■A[n\, 


where  we  assumed  that  g^,  gj  (and  ,  1^)  arc  linearly  independent,  and 


(4.51) 


(4.52) 


(4.53) 


(4.54) 


(4.55) 


/ 

T  i-n 

r  1 

J  T  1 

= 

T 

•[go  glj 

^  , 

\ 

Ls  J 

y 

Erom  (4.55)  the  two  information  symbols  are  given  by 

m[n]  =  }^  ■  ^n\, 


(4.56) 


(4.57) 
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and 


m\n-\\  =  1^  m[n]  =  •  A[«  +  l]. 

A  combination  of  (4.57)  and  (4.58)  yields 


g  ■  X[n\  =  g  ■  X[n  +  \'\<^g  ■  X[n\ - g  ■  X[n  + 1]  =  0, 


and  in  matrix  form 


[jg  -g~  ■ 


X[n\ 
A[n  +  1] 


=  0. 


(4.58) 


(4.59) 


(4.60) 


In  order  to  caleulate  k  = 


n 

obvious  that 


A[n] 
A[n  +  \] 


■  g[n] 


-k^ 


,  we  start  by  calculating  the  autoeorrelation  ma- 


^^[n  +  1]]  with  help  of  Equation  (4.53).  From  (4.60)  it  is 


fc"  -fe']'-RA.,=0.  (4.61) 

Therefore,  k  is  equal  to  the  eigenvector  corresponding  to  the  minimum  eigenvalue  of 
.  Once  the  estimate  for  k  is  known,  the  estimated  sequence  m\n\  is  derived  from 
(4.55).  The  equalization  algorithm  mentioned  above  is  tailored  for  the  speeific  ease  that, 
at  any  given  time,  y\n\  contains  only  two  symbols  (slight  misalignment).  If  this  is  not  the 

case  and  y\n\  eontains  more  than  two  symbols,  the  same  approach  can  be  followed  (in 

this  case  there  are  more  than  two  dominant  eigenvalues)  that  would  lead  to  a  similar  algo¬ 
rithm. 

C.  PROPOSED  DS-SS  SCHEME 

The  combination  of  the  DS-SS  modulation  seheme  and  the  receiver  structure 
based  on  the  aforementioned  blind  equalization  algorithms  was  proven  very  robust  as  part 
of  a  previous  thesis  [2].  However,  there  is  a  major  problem  associated  with  this  seheme. 
The  desired  subspaces  can  only  be  estimated  offline  since  the  whole  received  packet  is 
needed  in  order  to  calculate  the  data  correlation  matrix  7?^.  This  places  a  computational 

burden  that  increases  as  the  paeket  length  inereases. 
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In  this  thesis  the  signal  and  noise  subspaees  are  estimated  online  (reeursively) 
with  help  of  the  CRV  updating  algorithm  presented  in  Chapter  IV.  By  following  this  ap- 
proaeh,  the  desired  subspaees  are  effieiently  estimated  in  a  reeursive  manner.  This  is 
done  by  feeding  eaeh  eolumn  of  the  reeeived  data  matrix  y\n\  (4.40)  into  the  CRV  algo¬ 
rithm.  The  eolumns  of  y\n\  represent  the  reeeived  spread  bits,  eorresponding  to  one 
information  bit  spread  by  the  PN-eode  and  eorrupted  by  the  fading  ehannel. 

In  the  previously  proposed  eommunieation  seheme,  the  ehannel  impulse  response 
was  generated  using  the  Bellhop  underwater  aeoustie  propagation  model  [5].  The  time 
variability  of  the  ehannel  was  introdueed  to  the  statie  impulse  response  by  allowing  eaeh 
of  the  multipath  eomponents  to  vary  randomly  within  an  arbitrary  range.  The  random 
phase  shift  eaused  by  the  miero-paths  was  approximated  by  adding  a  random  phase  faetor 
to  eaeh  of  the  amplitudes. 

In  this  thesis  the  ehannel  impulse  response  is  modeled  using  the  Monterey-Miami 
parabolie  equation  model  (MMPE)  [7]  deseribed  in  Chapter  III.  The  broadband  impulse 
response  is  generated  for  the  frequeney  range  11, 500  ±1024  Hz  (BW  =  2048  Hz)  under 

various  perturbations  ineluding  the  influenee  of  Doppler  due  to  souree  motion.  The 
worst-ease  seenario  is  examined  by  evaluating  the  CRV  algorithm  using  a  time -variable 
impulse  response  that  also  aeeounts  for  the  interfaee  roughness,  the  turbulenee,  the  inter¬ 
nal  waves,  and  the  volume  seattering  eommonly  found  in  shallow  water  environments. 

D,  SIMULATION  RESULTS 

The  performanee  of  the  CRV  updating  algorithm  evaluated  under  different  ehan¬ 
nel  eonditions  and  design  parameters  (paeket  size,  bit  rate,  ete.)  is  examined  next.  The 
DS-SS  system  strueture  of  Figure  31  was  implemented  in  MATLAB  and  the  simulation 
results  are  presented  in  terms  of  Bit-Error-Rate  (BER)  vs.  Bit  Energy  per  Noise  Power 
(Ej/A^).  Eaeh  simulation  involved  sending  an  adequate  number  of  paekets  for  eaeh 
value  ofE^/Vg  in  order  to  ereate  suffieient  statisties.  The  is  a  sealed  version  of 

the  signal-to-noise  ratio  (SNR)  defined  as 

F  BW 

^  =  SNR-^,  (4.62) 

N,  R, 
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where  is  the  available  communication  bandwidth  and  is  the  bit  rate.  The  range  of 

the  desired  used  for  the  various  simulations  was  between  0  and  30  dB.  The  BER 

(probability  of  a  bit  error)  is  a  metric  used  to  describe  the  reliability  of  a  digital  commu¬ 
nication  system  quantitatively  and  is  defined  as 

Total  Number  of  Errors 

BER= - .  (4.63) 

Total  Number  of  Bits  Sent 

1,  Performance  in  AWGN 

The  received  sequence  is  s\n\  +  M\n\  where  s\n\  is  the  spread  information  se¬ 
quence  andw[n]  is  the  AWGN.  In  this  case  the  channel  corruption  is  not  taken  into  ac¬ 
count;  thus,  the  results  of  this  simulation  can  be  treated  as  a  benchmark  for  evaluating  the 
receiver  performance.  Eigure  34  shows  the  resulting  BER  for  different  packet  sizes  and 
Eigure  35  the  BER  for  various  data  rates  . 


Eigure  34.  BER  in  an  AWGN  Channel  When  Varying  the  Packet  Size 
{R,,=  40bps,  =  2400  cps) . 
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Figure  35.  BER  in  an  AWGN  Channel  When  Varying  the  Data  RateR^ 

(packet  size  =  72 bits,  =2400cps) . 

From  the  results  presented  above  two  observations  can  be  made.  First,  we  see 
that,  if  longer  transmission  packets  are  used,  the  performance  is  not  degraded  in  this  ideal 
case.  The  same  conclusion  stands  for  the  case  where  a  higher  data  rate  is  used.  The  per¬ 
formance  seems  to  be  better  for  a  packet  of  144  bits  (Figure  34)  and  for  a  data  rate  of  60 
bps  (Figure  35).  This  is  possibly  because  of  statistical  variations. 

2,  Performance  in  Additive  Color  Noise 

The  colored  noise  is  produced  by  filtering  white  noise  by  the  7*^-order  FIR  filter 
given  by 


0 


Qiz)  = 


1 


z*’  -0.9z'  +0.4z"  -0.2z'  +0.1z"  -O.lz  ' 

This  is  the  same  filter  used  in  [2]  so  that  the  performance  results  may  be  compared. 


(4.64) 


Figure  36  illustrates  the  performance  of  the  new  subspace  decomposition  algo¬ 
rithm  (CRV)  in  additive  colored  noise  when  varying  the  packet  size.  The  results  of  this 
simulation  are  worse  than  the  white  noise  case  by  approximately  4  dB.  Figure  37  demon¬ 
strates  the  receiver  performance  for  a  72-bit  packet  when  the  bit  rate  is  varied  in  the  range 
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of  40  to  100  bps.  It  is  shown  once  again  that  the  algorithm  is  stable  even  for  a  bit  rate  of 
100  bps. 


Figure  36.  Performanee  in  Additive  Color  Noise  When  Varying  the  Paeket  Size 

{Ri,=  40bps,  =  2400  eps) . 


Figure  37.  Performanee  in  Additive  Color  Noise  When  Varying  the  Data  Rate 
(paeket  size  =  72 bits,  R^  =2400  eps) . 
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3,  Performance  in  Multipath  and  AWGN 

The  received  sequence  is  s{n\^h{n\  +  M\n\  where  s{n\  is  the  spread  by  the 
PN-code  data  sequence,  >v[«]  is  the  AWGN,  and  h{n\  is  the  channel  impulse  response. 
The  impulse  responses  used  for  the  simulations  are  the  ones  computed  with  help  of  the 
MMPE  model.  The  performance  of  the  CRV  algorithm  was  evaluated  for  all  the  different 
perturbations  described  in  Chapter  III  and  the  results  are  shown  in  Figure  38  through  Fig¬ 
ure  44. 

As  already  mentioned,  the  communication  scheme  of  Figure  3 1  is  tailored  for  a 
channel  impulse  response  equal  to  the  processing  gain  N.  For  all  the  following  simula¬ 
tions,  =  40 bps  and  =2400cps ,  thus  N  is  60.  This  is  the  reason  all  results  are  pre¬ 
sented  for  a  channel  length  of  at  least  60  time  samples. 

Among  the  individual  perturbations,  the  most  challenging  are  those  corresponding 
to  the  interface  roughness  (Figure  39)  and  turbulence  (Figure  42).  In  those  cases  the  re¬ 
sults  are  good  only  for  a  channel  length  of  60  time  samples. 

When  all  the  examined  perturbations  are  combined  (Figure  43)  the  results  are 
poor  and  on  the  average  worse  than  the  previous  cases  by  1 8  dB  (for  an  impulse  response 
of  60  samples).  The  results  of  the  worst-case  scenario  where  the  Doppler  is  also  taken 
into  account  are  presented  in  Figure  44.  In  this  case,  equalization  cannot  be  performed, 
even  for  a  channel  length  of  60  time  samples.  However,  it  must  be  noted  that  the  exam¬ 
ined  perturbations  are  expected  to  be  realistic  but  strong  perturbations.  It  is  doubtful  that 
all  of  them  combined  will  be  encountered  in  many  shallow  water  environments. 
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Figure  38.  No  Perturbation;  BER  in  a  Fading  Channel  with  AWGN  When  Varying 
the  Channel  Length  (paeketsize  =  72 bits,  =40 bps,  =2400eps) . 


Figure  39.  Interface  Roughness  Perturbation;  BER  in  a  Fading  Channel  with  AWGN 
When  Varying  the  Channel  Length  (paeketsize  =  72 bits, =40 bps,  =2400cps) . 
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Figure  40.  Internal  Wave  Perturbation;  BER  in  a  Fading  Channel  with  AWGN  When 
Varying  the  Channel  Length  (packet  size  =  72  bits,  =40  bps,  =2400cps) . 


Figure  41 .  Bottom  Volume  Perturbation;  BER  in  a  Fading  Channel  with  AWGN 
When  Varying  the  Channel  Length  (packet  size  =  72  bits,  =40  bps,  =2400cps) . 
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Figure  42.  Turbulence  Perturbation;  BER  in  a  Fading  Channel  with  AWGN  When 
Varying  the  Channel  Length  (packet  size  =  72  bits,  =40 bps,  =2400cps) . 


Figure  43.  All  Perturbations  Combined  (no  Doppler):  BER  in  a  Fading  Channel  with 
AWGN  When  Varying  the  Channel  Length 
(packet  size  =  72 bits,  =40 bps,  =2400  cps) . 
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Figure  44.  All  Perturbations  Combined  (with  Doppler):  BER  in  a  Fading  Channel 
with  AWGN  When  Varying  the  Channel  Length 
(packet  size  =  72 bits,  =40 bps,  =2400  cps) . 

4,  Performance  in  Multipath  and  Additive  Color  Noise 

This  is  the  same  case  as  the  previous  one  but  colored  noise  was  added  instead  of 
white.  The  equalization  algorithm  used  in  this  case  is  described  by  Equations  (5.11) 
through  (5.15)  where  once  again  the  CRV  is  used  to  estimate  the  signal  and  noise  sub¬ 
spaces.  The  resulting  BER  for  a  channel  length  of  60  time  samples  is  shown  in  Figure  45. 
The  algorithm  works  well  for  all  the  individual  perturbations  but  not  when  all  of  them  are 
combined. 


77 


Figure  45.  Performance  in  a  Fading  Channel  with  Additive  Color  Noise  When  Vary¬ 
ing  the  Channel  Length  (packet  size  =  72  bits,  =40 bps,  =2400cps) . 

5,  Performance  When  Synchronization  is  Lost 

When  the  PN-code  in  the  transmitter  and  its  replica  in  the  receiver  are  synchro¬ 
nized  (aligned  in  time),  the  received  sequence  is  despread  correctly.  However  due  to  the 
volatile  underwater  environment,  this  may  not  always  be  the  case.  When  synchronization 
is  lost,  the  dominant  eigenvalues  of  the  received  covariance  matrix  are  more  than  one. 
When  at  any  given  time  y[n]  contains  two  symbols  (the  current  symbol  m[n]  and  a  por¬ 
tion  of  the  previous  symbol  m[n-\]).  Equations  (4.49)  through  (4.61)  may  be  used  for 
the  recovery  of  the  transmitted  sequence.  Therefore,  in  order  to  simulate  the  loss  of  syn¬ 
chronization,  a  random  delay  (of  duration  up  to  half  a  bit)  was  added  to  the  received  sig¬ 
nal.  The  performance  results  of  the  algorithm  are  poor  for  all  individual  perturbations. 
For  the  case  in  which  all  perturbations  are  combined  and  subsequently  the  Doppler  Effect 
is  added  the  resulting  BER’s  are  shown  in  Eigure  46. 
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Figure  46.  BER  for  the  Case  in  whieh  Synchronization  of  the  PN-code  is  Lost. 

6,  Summary  of  Performance 

In  this  chapter,  the  performance  of  the  three  different  equalization  algorithms, 
with  the  CRV  used  for  the  estimation  of  the  signal  and  noise  subspace,  was  examined  un¬ 
der  various  design  parameters  and  channel  conditions. 

When  synchronization  between  the  spreading  sequence  and  its  replica  is  achieved, 
the  results  of  the  proposed  communication  scheme  are  good  for  all  the  individual 
perturbations.  For  a  fading  channel  with  AWGN,  a  BER  of  10'^  is  possible  for 

>  13dB  whereas  in  a  fading  channel  with  additive  color  noise  the  same  BER  is 

achieved  forE^/A^  >  17dB  . 

Among  the  examined  perturbations,  the  ones  corresponding  to  the  interface 
roughness  and  the  turbulence  have  more  impact  on  the  performance  of  the  communica¬ 
tion  scheme  than  the  others.  The  small-scale  fluctuations  of  sound  speed  within  the  water 
column  caused  by  turbulence,  and  the  scattering  of  the  acoustic  energy  interacting  with 
the  rough  bottom  create  a  very  challenging  environment  for  underwater  wireless  commu¬ 
nications  due  to  the  greater  breakdown  in  the  coherence  of  the  propagating  field. 
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The  performance  of  the  algorithms  is  not  the  desired  one  when  all  the  perturba¬ 
tions  are  combined  and  Doppler  is  added  to  account  for  the  source  motion.  But  as  already 
mentioned  the  modeled  perturbations  are  rather  strong,  creating  a  shallow  water  envi¬ 
ronment  more  challenging  than  usual. 

When  synchronization  is  not  achieved  and  the  algorithm  described  earlier  is  used 
for  the  recovery  of  the  information  sequence,  the  results  are  very  poor.  As  shown  in 
Figure  46,  a  BER  <  10'^  cannot  be  achieved  even  forE^/A^  =  25dB  .  In  the  last  chapter, 
the  important  findings  of  this  thesis  research  are  summarized  and  possible  areas  for  fol- 
low-on  work  are  presented. 
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VI.  CONCLUSION 


The  primary  purpose  of  this  thesis  was  to  evaluate  the  performance  of  the  CRV 
updating  algorithm  developed  in  a  previously  proposed  communication  scheme  intended 
for  use  in  an  underwater  acoustic  network  operating  in  a  shallow-water  environment.  The 
transmitter  and  receiver  were  implemented  in  MATLAB  and  the  underwater  channel  was 
modeled  using  the  Monterey-Miami  Parabolic  Equation  model  (MMPE).  Various  simula¬ 
tions  were  performed  under  different  channel  conditions  (interface  roughness,  turbulence, 
etc.)  and  design  parameters  (packet  size,  bit  rate). 

A,  FINDINGS 

The  CRV  was  proposed  [9]  as  an  alternative  to  the  eigenvalue  decomposition 
(EVD).  The  main  advantage  offered  by  the  CRV  is  that  it  can  be  updated  online  at  less 
computational  cost  compared  to  the  EVD.  As  shown  in  Chapter  IV,  the  estimates  of  the 
signal  and  noise  subspace  produced  by  the  recursively  updated  algorithm  are  very  accu¬ 
rate. 

In  order  to  examine  the  performance  of  the  CRV,  three  distinct  cases  were  simu¬ 
lated;  a  shallow-water  fading  channel  with  AWGN,  a  fading  channel  with  .colored  noise, 
and  the  case  in  which  synchronization  between  the  receiver  and  the  transmitter  were  not 
achieved.  Under  all  individual  perturbations  described  in  Chapter  III,  the  equalization 
algorithm  performed  well  and  achieved  a  BER  of  <10^^ for  values  of  >  13  dB 

andEj/V^  >  17dB  for  white  and  colored  noise,  respectively. 

Among  the  modeled  perturbations,  the  ones  corresponding  to  the  interface  rough¬ 
ness  and  turbulence  were  the  most  challenging  for  the  proposed  equalization  algorithm. 
As  mentioned  in  Chapter  III,  these  perturbations  cause  the  greater  breakdown  in  the  co¬ 
herence  of  the  propagating  field  resulting  in  a  severe  multipath  environment. 

On  the  other  hand,  when  all  perturbations  were  combined  and  the  Doppler  effect 
due  to  source  motion  was  considered,  the  resulting  bit-error-rates  were  not  acceptable. 
Eor  the  case  of  no  synchronization  between  the  spreading  and  de-spreading  sequence,  the 
performance  of  the  acquisition  algorithm  was  also  very  poor. 
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B,  FUTURE  WORK 

The  channel  impulse  response  was  obtained  by  post  processing  the  output  of  the 
MMPE  acoustic  propagation  model,  which  represented  the  response  of  the  simulated  un¬ 
dersea  environment.  Therefore,  the  proposed  communication  scheme  was  tested  under 
simulated  conditions.  It  would  be  very  interesting  to  evaluate  its  performance  using  real 
impulse  responses  obtained  from  experimental  data  in  a  variety  of  shallow-water  chan¬ 
nels. 

In  the  current  implementation,  the  desired  subspaces  (signal,  noise)  were  esti¬ 
mated  online  with  help  of  the  recursively  updated  CRV  algorithm.  However,  the  equali¬ 
zation  based  on  matched  filter  theory  is  still  performed  offline.  More  research  on  design¬ 
ing  an  effective  communication  scheme  that  could  operate  online  is  needed.  This  can  be 
partially  achieved  by  using  training  bits  in  each  transmission  or  by  using  the  flrst  packet 
as  a  training  transmission. 

The  CRV  updating  algorithm  could  be  even  more  computationally  efficient  if  the 
Lanczos  bidiagonalization  algorithm  were  used  to  compute  the  eigenvalues  of  the'P  ma¬ 
trix  given  by  Equation  (4.18). 

There  is  a  family  of  rank-revealing  orthogonal  decompositions,  referred  to  as 
UTV  decompositions  [33],  offered  as  alternatives  to  the  SVD.  These  algorithms  are  able 
to  provide  all  the  necessary  subspace  information,  and  most  importantly,  can  be  updated 
online.  A  comparison  of  their  performance,  with  the  results  obtained  with  the  CRV  de¬ 
composition,  could  lead  to  the  design  of  a  better  signal  processing  algorithm  for  the  re¬ 
covery  of  the  information  signal  in  a  multipath  channel. 

Einally,  the  performance  of  the  proposed  transmitter/receiver  structure  should  be 
evaluated  as  part  of  a  complete  communication  system  with  error  correction  coding  and 
modulation. 
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APPENDIX. 


MATLAB  CODE 


The  MATLAB  source  code  used  for  the  simulation  of  the  communication  system 
structure  of  Figure  31  and  the  implementation  of  the  CRV  updating  algorithm  is  pre¬ 
sented  in  this  appendix.  The  overall  controlling  m-file  is  called  CRV_simulations.m.  All 
the  other  necessary  functions  (m-files)  are  called  by  running  this  file  and  are  listed  bel¬ 
low. 

CRV.m  -  This  function  implements  the  CRV  updating  algorithm  described  in 
Chapter  IV.  The  output  of  this  function  is  the  recursively  updated  correlation  matrix  0  , 
and  the  matrices  W  and  V  needed  for  the  initialization  of  the  next  step  of  the  algorithm. 

up.m  -  This  function  is  used  to  upsample  the  message  signal  to  the  chipping  rate 
so  that  the  PN-code  may  be  appended. 

spredesp.m  -  This  function  is  used  to  spread  the  binary  message  signal  with  the 
Gold  code,  or  de-spread  the  received  sequence  with  the  replica  of  the  Gold  code  used  at 
the  transmitter. 

power _of_noise.m  -  This  function  calculates  the  noise  power  associated  with  a 
given  [dB],  a  bit  length  [s],  and  a  sampling  period  [Hz].  The  output  is  used 

to  generate  the  white  or  colored  noise  corresponding  to  the  noise  power. 

house. m  -  This  function  computes  the  Householder  matrix  P  required  to  clear 
specific  values  of  the  received  sequence  (a  required  step  in  the  CRV  updating  algorithm). 
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%  Overall  Controlling  Program 

%  T/R  SIMULATION 

%  Note:  all  seeondary  variables  and  needed 
%  funetions  are  ealled  by  exeeuting  this  m-file 

% 

%  developed  by  Pavlos  Angelopoulos  Mar  2004 
%  last  modified  4/13 

elear; 

ele; 

%  Offer  options  for  different  simulations 
disp('Choose  one  of  the  following  options:'); 
disp('l.  White  Noise  only'); 
disp('2.  Color  Noise  only'); 
disp('3.  Fading  Channel+AWGN'); 
disp('4.  Fading  Channel+ACN'); 
disp('5.  Loss  of  Reeeiver  Synehronization'); 

option=input(' '); 
dispC  '); 

if  option==3|option==4; 

%  Offer  options  for  different  perturbations 

disp('Choose  one  of  the  following  perturbations:'); 

disp('L  Basie'); 

disp('2.  Interfaee  Roughness'); 

disp('3.  Internal  Wave'); 

disp('4.  Volume'); 

disp('5.  Turbulenee'); 

disp('6.  All  Perturbations  (without  Doppler)'); 
disp('7.  All  Perturbations  (with  Doppler)'); 
pert=input(' '); 
dispC '); 

%  Load  the  desired  Impulse  Response 
if  pert==l;  load  ImpRes  basie; 
elseifpert==2;  load  ImpRes  intrough; 
elseifpert==3;  load  ImpRes  intwave; 
elseifpert==4;  load  ImpRes  volume; 
elseifpert==5;  load  ImpRes  turbulenee; 
elseifpert==6;  load  ImpRes  allpert; 
elseifpert==7;  load  ImpRes  doppler; 
end  %  end  of  "if  pert" 
dispC '); 

end  %  end  of  "if  option==3|option==4" 

if  option==3|option==4|option==5; 
disp('Choose  the  Channel  Length:'); 
ehannel_length=input(' '); 
dispC '); 

end  %  end  of  "if  option==3|option==4|option==5" 


disp('Choose  #  of  Simulations:'); 
trials=input(' '); 
dispC  '); 
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%■ 


70 . 

Rb  =  40; 

%  bit  rate  [bps] 

Re  =  2400; 

%  chip  rate  [cps] 

BW=2*Rc; 

%  bandwidth  of  spread  signal 

Tb  =  1/Rb; 

%  bit  Period  [s] 

Tc  =  1/Rc; 

%  chip  Period  [s] 

N  =  Rc/Rb; 

%  processing  gain 

P=Tb/Tc; 

%  Upsampling  required  to  match  chip  sequence 

load  pncode; 

%  load  PN  chipping  sequence  "pn  code.mat" 

packet  length  =  72; 

%  number  of  bits  transmitted  per  packet 

number  of  packets=l; 

%  number  of  packets  transmitted 

EbNo  dB=hnspace(l,25,12);  %  determine  Eb/No  values  to  run  simulations 

disp('START'); 

dispC '); 

% - 

% - 

% 

TRANSMITTER 

% - 

% - 

for  simulation=l:trials; 

%  #of  simulations 

txt=sprintf('Simulation  #%3  d',  simulation) ; 
disp(txt); 

tic; 

%  starts  clock  to  measure  simulation  run  time 

txt2=sprintf('Simulation  #%3d', simulation); 
h=waitbar(0,txt2); 

for  11=1  :length(EbNo_dB);  %  run  simulation  for  each  of  the  EbNo  values 
waitbar(n/length(EbNo_dB),h); 


% - 

%  Generate  the  binary  data  sequence  of  +1,-1 

% - 

data_seq=sign(randn(l,number_of_packets*packet_length)); 


% - 

%  Upsample  by  P=Tb/Tc  to  match  the  chipping  rate 

% - 

upsampled_data_seq=up(data_seq,P); 

m=length(upsampled_data_seq); 


% - 

%  Spread  the  data  using  with  chipping  sequence  c[n] 

% - 

spreaded_signal=conv(upsampled_data_seq,pn_code(l:P)); 
spreaded_signal=spreaded_signal(  1  :m) ; 


% - 

%  Calculate  the  noise  power  associated 
%  with  the  simulated  EbNo  values 
% - 

noise_power=power_of_noise(spreaded_signal,Tb,l/BW,EbNo_dB(n)); 
noise  std  =  sqrt(noise  jiower);  %  noise  standard  deviation 
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% - 

% - 

%  CHANNEL 

% - 

% - 


ifoptioii==l; 

white_noise=randn(  1  ,length(spreaded_signal)) ; 

received_seq=  spreaded  signal  +  white_noise*noise_std;%  Add  White  Noise 
elseif  optioii==2; 
white_noise=randn(  1  ,m) ; 

color_noise=filter(l,[l  -0.9  0.4  -0.2  0.1  -0.1],white_noise);%  fdter  white  noise  to  produce  color  noise 
received_seq=  spreaded  signal  -i-  color_noise*noise_std;%  Add  Color  Noise 
elseif  option==3; 

spreaded_signal=conv(imp_res(l:channel_length),spreaded_signal); 
white_noise=randn(  1  ,length(spreaded_signal)) ; 

received_seq=  spreaded  signal  -i-  white_noise*noise_std;%  Add  White  Noise 
elseif  option==4; 

spreaded_signal=conv(imp_res(l:channel_length),spreaded_signal); 
white_noise=randn(  1  ,length(spreaded_signal)) ; 
color_noise=fdter(l,[l  -0.9  0.4  -0.2  0.1  -0.1],white_noise); 
received_seq=  spreaded  signal  -i-  color_noise*noise_std;%  Add  White  Noise 
else 

load  ImpResdoppler; 

spreaded_signal=conv(imp_res(l:channel_length),spreaded_signal); 
white_noise=randn(  1  ,length(spreaded_signal)) ; 

received_seq=  spreaded  signal  -i-  white_noise*noise_std;%  Add  White  Noise 
end  %  end  of  "if  option==l " 


% - 

% - 

%  RECEIVER 

% - 

% - 


if  option==5; 

%  Add  random  time  tO  uniformly  distributed  between  0  and  P 
%  to  simulate  the  Loss  of  Synchronization  between  the  codes 
tO=round(rand*P/4)-l-l ; 

received_seq=received_seq(tO:length(received_seq)); 

end 

% - 

%  Despread  the  received  sequence  using  a  replica  of  c[n] 

% - 

despreaded_signal=spredesp(received_seq(l  :m),pn_code(l  :P)); 
despreaded_signal=despreaded_signal(l:m); 


% - 

%  Delay  by  z'^-l  and  Downsample  by  P 

% - 

Y =re  shape  (de  spreadedsignal,  P,packet_length) ; 
columns=size(Y,2); 
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% - 

%  Implement  the  CRV  algorithm 
% - 


%  Initialization 
W=zeros(P,P); 

V=eye(P); 

1^0;  %  Rank  of  matrix  W  (r=0) 
fork=l:columns; 

[phi_crv,W,V,r]=CRV(Y(:,k),W,V,r);  %  call  the  "CRV.m"  function 
end 


% - 

%  Compute  the  matched  filter  coefficients 
% - 


if  option==  1  |option==3 ; 

[EigVec,EigVal]  =  eig(phi_crv); 

Opt_filter=EigVec(:,l);  %  Eigenvector  of  the  largest  eigenvalue 

elseif  option==2|option==4; 

corrmatrix_color_noise=color_noise*color_noise'; 
[EigVec_color,EigVal_color]=eig(corrmatrix_color_noise); 
R=inv(sqrtm(EigV  alcolor)) ; 

R_sq_inv=EigVec_color*R*EigVec_color'; 

R_mah=R_sq_inv  *phi_crv  *  R_sci_inv ; 

[Eig  V  ec_mah,Eig  Val_mah]=eig(R_mah) ; 

Opt_filter=R_sq_inv  *  Eig  V  ec_mah( : ,  1 ) ; 
elseif  option==5; 

[EigVec,EigVal]  =  eig(phi_crv); 

Opt_filter=EigVec(:,l);  %  Eigenvector  of  the  largest  eigenvalue 

Eig  Val=diag(EigV  al) ; 

if  EigVal(2)>0.2*EigVal(l);  %  Assume  there  is  no  synchronization 
lamda=EigVec(:,  1 :2)'*  Y ; 

E=[lamda( : ,  1  :length(lamda)- 1 ) ;  lamda( :  ,2 :  length(lamda))] ; 

R_lamda=E*E'; 

[Ve,  De]=eig(R_lamda); 
lam_e=diag(De); 

I=find(lam_e==min(abs(lam_e))); 

ke=Ve(:,I); 

kl=ke(l:2); 

k2=-ke(3:4); 

end 

end 


% - 

%  Filter  the  received  signal 

% - 

if  option— 5  &  EigVal(2)>0.1*EigVal(l); 

decoded_seq=sign([k  1  ',k2']  *E) ; 
else 

decoded_seq=sign(Opt_filter'*Y) ; 
end; 
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% - 

%  Calculate  Errors 
% - 

if  option— 5  &  EigVal(2)>0.1*EigVal(l); 
diff=length(decoded_seq)-length(data_seq); 
ifdifE>0 

err=length(find(data_seq-decoded_seq(diff:length(data_seq)+diff-l)~=0)); 
elseif  diff<0 
diff=abs(diff); 

err=length(find(data_seq(diff:lengtli(decoded_seq)+diff-l)-decoded_seq~=0)); 

else 

err=sum(abs(data_seq-decoded_seq))/2;  %  #  of  errors 
end  %  end  of  "if  diff>0" 
else 

err=sum(abs(data_seq-decoded_seq))/2;  %  #  of  errors 
end  %  end  of  "if  option— 5  &  EigVal(2)>0. 1  *EigVal(l)" 
if  err>packet_length/2; 

channel_errors(n)=packet_length-err; 

else 

channel_errors(n)=err; 
end  %  end  of  "if  err>packet_length/2" 

end  %  end  of  n  loop 

close(h); 

toe; 

summary(simulation,:)=channel_errors./packet_length; 

end  %  end  of  simulations  loop 

avg_errors=mean(summary) ; 

result=[EbNo_dB;  avg  errors]; 

% - 

%  Plot  Simulation  Results 

% - 

result(2,find(result(2,:)— 0))=10^-6; 
figure; 

semilogy(result(l,:),result(2,:),'-VMarkerSize',8); 
grid  on; 

ylim([10^-5  l]);xlim([l  30]); 
ylabel('BERVfontsize',  1 3); 
xlabel('Eb/No  (dB)Vfontsize',13); 
legend(sprintf('Packet  size:  %3.0fbits',packet_length)); 

disp('STOP'); 

%  offer  option  to  save  the  simulation  results 
%  for  post  processing  (plots) 

savedata=input('Save  Simulation  Results?  (y  or  n)  Vs'); 
dispC '); 

if  savedata  —  'y'  |  savedata  ==  'Y' 
fileout=input('Enter  Output  filename  (without  extension):  ','s'); 
eval(['save  '  fileout '  result  trials  packet  length  Re;']); 
end 

if  savedata  —  'n'  |  savedata  ==  'N' 
end 
break 

o/„ - end  of  Code - 
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function  [phi_crv,W,V,r]=CRV(signal,W,V,r) 

%  CRV  UPDATED  ALGORITHM 

%  "signal"  is  the  new  signal  used  to  update  the  correlation  matrix 
%  "  W"  and  "V"  are  the  matrices  needed  to  initialize  the  algorithm 
%  "r"  is  the  rank  of  W 

%  "phi  crv"  is  the  esimate  of  the  correlation  matrix 

% 

%  developed  by  Pavlos  Angelopoulos  Feb  2004 
%  last  modified  3/28 

a=l;  %  Fading  Factor 
P=length(signal); 

W=a*W; 
z=V'*  signal; 


ifr~=0; 
ifr— P; 

W=W+z*z'; 

else 

P_hat=house(z(r+l  :end));  %  Householder  vector:  keep  r+1  element  and  zero  the  rest 

P  diag=[eye(r)  zeros(r,P-r);zeros(P-r,r)  P  hat];  %  P=diag(Ir,P  hat) 


W=P_diag'*(W+z*z')*P_diag; 

V=V*P_diag; 

warning  off; 

eigenv=abs(eigs(  W  ( 1  :r+ 1 , 1  :r+ 1  ),r+ 1 )) ; 
if  eigenv(r+l)>0.1*eigenv(r); 
r=r+l; 
end 
end 

phi_crv=V*W*V'; 

else 

W=W+z*z'; 

phi_crv=V*W*V'; 

r=rank(W); 

end 

o/„ - end  of  function  CRV - 


%  update  matrix  W 
%  update  matrix  V 

%  compute  eigenvalues  of  (r+l)x(r+l)  block  of  W 
%  Threshold  to  check  rank  change 


function  y=up(x,M) 

%  Upsample:  A  function  to  upsample  a  signal  x  by  an  integer  M 

% 

%  developed  by  Pavlos  Angelopoulos  Nov  2003 
%  last  modified  3/28 

n=l:length(x); 

y(n*M)=x; 

y=[y(M:end)  zeros(l,M-l)]; 

% - End  of  function  "up" - 
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function  noise_power=power_of_noise(x,Tb,Ts,EbNo_dB) 

%  Noise  Power 

%  A  function  to  calculate  the  noise  power  associated  with  a 
%  given  Eb/No  (dB),  bit  length  Tb,  and  sampling  period  Ts 

% 

%  developed  by  Pavlos  Angelopoulos  Nov  2003 
%  last  modified  1 1/03 

signal_power=sum(x.^2)/length(x); 

EbNo=10^(EbNo_dB/10); 

Eb=signal_power*Tb;  %  energy  per  bit 

No=Eb/EbNo; 

sigma_No=sqrt(No/(2*Ts));  %  std  deviation  of  noise 
noise_power=sigma_No''2; 

% - End  of  function  "power  of  noise" - 


function  P=house(x) 

%  A  function  to  compute  the  householder  matrix  P  such  that 
%  transformed_x=P*x  is  a  vector  whose  elements  except  the 
%  first  one  are  equal  to  0. 

%  X  must  be  defined  as  a  column  vector 
% 

%  developed  by  P.Angelopoulos  Nov  2003 
%  last  modified  1 1/24 

k=length(x); 
el=[l,zeros(l,k-l)]'; 
v=x+sign(x(  1  ))*norm(x,2)*e  1 ; 

P=eye(k)-2*v*v'./(v'*v); 

% - End  of  function  "house" - 


function  data=spredesp(x,code) 

%  A  function  to  spread  (chip)  a  binary  signal  x  with  gold  code  "code" 
%  or  de-spread  the  received  signal  by  a  replica  of  the  "code" 

% 

%  developed  by  P.Angelopoulos  Oct  2003 
%  last  modified  10/23 

k=length(x); 

m=length(code); 

repeated_code=repmat(code,  1  ,ceil(k/m)); 
data=x.  *repeated_code; 

% - End  of  function  "spredesp" - 


90 


LIST  OF  REFERENCES 


[1]  P.  Duke,  “Direct-Sequence  Spread-Spectrum  Modulation  for  Utility  Packet 
Transmission  in  Underwater  Acoustic  Communication  Networks”,  Master’s  The¬ 
sis,  Naval  Postgraduate  School,  Monterey,  California,  2002. 

[2]  G.  N.  Pelekanos,  “Performance  of  Acoustic  Spread-Spectrum  Signaling  in  Simu¬ 
lated  Ocean  Channels”,  Master’s  Thesis,  Naval  Postgraduate  School,  Monterey, 
California,  2003. 

[3]  J.  A.  Rice,  R.  K.  Creber,  C.  L.  Fletcher,  P  .A.  Baxley,  D  .E.  Rogers,  and  D.  C. 
Davsion,  “Seaweb  Underwater  Acoustic  Nets,”  Space  and  Naval  Warfare  Systems 
Center  San  Diego  Biennial  Review  2001,  pp.  234-243,  2001. 

[4]  Joe  Rice,  Bob  Creber,  Chris  Fletcher,  Paul  Baxley,  Ken  Rogers,  Keyko  McDon¬ 
ald,  Dave  Rees,  Michael  Wolf,  Steve  Merriam,  Rami  Mechio,  John  Proakis,  Ken 
Scussel,  Dave  Porta,  John  Baker,  Jim  Hardiman,  and  Dale  Green,  “Evolution  of 
Seaweb  Underwater  Acoustic  Networking”,  IEEE,  pp.  2007-2017,  2000. 

[5]  M.  B.  Porter,  “Ocean  Acoustics  Eibrary,”  [http://oalib.saic.com],  last  accessed 
February  2004. 

[6]  NATO  SACEANT,  “Acoustic  Models”,  [http://www.saclantc.nato.int/frameset- 
ch6.html],  last  accessed  February  2004. 

[7]  K.  B.  Smith,  “Convergence,  stability,  and  variability  of  shallow  water  acoustic 
predictions  using  a  split-step  Fourier  parabolic  equation  model”,  J.  Comp. 
Acoust.,  9,  pp.  243-285,  2000. 

[8]  J.  G.  Proakis,  Digital  Communications,  4*  ed.,  pp.  126-121 ,  McGraw-Hill,  New 
York,  2001. 


[9]  E.  S.  Baker  and  R.  D.  DeGroat,  “A  Correlation-Based  Subspace  Tracking  Algo¬ 
rithm”,  IEEE  Transactions  on  Signal  Processing,  Vol.  46,  pp.  3112-3116,  1998. 


91 


[10]  Lawrence  E.  Kinsler,  Austin  R.  Frey,  Alan  B.  Coppers,  and  James  V.  Sanders, 
“Underwater  Acoustics,”  m  Fundamentals  of  Acoustics,  4*  ed.,  pp.  435-470,  John 
Wiley  &  Sons,  New  York,  2000. 

[11]  Finn  B.  Jensen,  William  A.  Kuperman,  Michael  B.  Porter,  and  Henrik  Schmidt, 
Computational  Ocean  Acoustics,  AlP  Press,  New  York,  1994. 

[12]  Robert  J.  Urick,  Principles  of  Underwater  Sound,  ed.,  McGraw-Hill,  New 
York,  1983. 

[13]  F.  M.  Brekhovskikh,  and  Yu.  P.  Lysanov,  Fundamentals  of  Ocean  Acoustics,  3'^‘^ 
ed..  Springer- Verlag,  New  York,  2002. 

[14]  Robert  J.  Urick,  Ambient  Noise  in  the  Sea,  Peninsula  Publishing,  Los  Altos,  Cali¬ 
fornia,  1986. 

[15]  Gordon  M.  Wenz,  “Acoustic  Ambient  Noise  in  the  Ocean;  Spectra  and  Sources,” 
Journal  of  the  Acoustical  Society  of  America,  Vol.  34,  No.  12,  pp.  1936-1956, 
1962. 


[16]  Rodney  F.  W.  Coates,  “Noise  and  Reverberation,”  in  Underwater  Acoustic  Sys¬ 
tems,  pp.  90-94,  John  Wiley  &  Sons,  New  York,  1989. 

[17]  www.benthos.com/pdf/Telesonar%20Transducers.pdf,  last  accessed  March  2004. 

[18]  Theodore  S.  Rappaport,  Wireless  Communications:  Principles  and  Practices,  2“‘^ 
ed.,  Prentice  Hall,  Upper  Saddle  River,  New  Jersey,  2002. 

[19]  F.  D.  Tapper!,  “The  parabolic  approximation  method”  (Chapter  V).  Lecture  Notes 
in  Physics,  Vol.  70,  Wave  Propagation  and  Underwater  Acoustics,  eds.  J.  B.  Kel¬ 
ler  and  J.  S.  Papadakis,  Springier- Verlag,  New  York,  1977. 

[20]  R.  H.  Hardin  and  F.  D.  Tapper!,  “Applications  of  the  split-step  Fourier  method  to 
the  numerical  solution  of  nonlinear  and  variable  coefficient  wave  equations,” 
SIAM  Rev.,  Vol.  15,  pg.  423,  1973. 


92 


[21]  J.  A.  Goff  and  T.  H.  Jordan,  “Stochastic  modeling  of  seafloor  morphology:  Inver¬ 
sion  os  Sea  Beam  data  for  seeond-order  statisties,”  J.  Geophys.  Res.,  Vol.  93,  pp. 
13589-13609,  1988. 

[22]  Robert  M.  Hill,  “Model-Data  Comparison  of  Shallow  Water  Aeoustie  Reverbera¬ 
tion  in  the  East  China  Sea,”  Master’s  Thesis,  Naval  Postgraduate  Sehool,  Mon¬ 
terey,  California,  2003 . 

[23]  T.  Yamamoto,  “Veloeity  variabilities  and  other  physieal  properties  of  marine 
sediments  measured  by  eroswell  aeoustie  tomography,”  Journal  of  the  Acoustical 
Society  of  America,  Vol.  98,  pp.  2235-2248,  1995. 

[24]  V.  I.  Tatarski,  Wave  Propagation  in  a  Turbulent  Medium,  Dover  Publications, 
London,  1961. 

[25]  T.  F.  Duda,  S.  M.  Flatte’,  and  D.  B.  Creamer,  “Modelling  Meter-Seale  Aeoustie 
Intensity  Fluctuations  From  Oceanic  Fine  Strueture  and  Microstructure,”  Journal 
of  Geophysical  Research,  Vol.  93,  No.  C5,  pp.  5130-5142,  1988. 

[26]  F.  S.  Henyey,  D.  Rouseff,  J.  M.  Groehoeinski,  S.  A.  Reynolds,  K.  L.  Williams, 
and  T.  E.  Ewart,  “Effects  of  Internal  Waves  and  Turbulenee  on  a  Horizontal  Ap¬ 
erture  Sonar,”  IEEE  Journal  of  Oceanic  Engineering,  Vol.  22,  No.  2,  pp.  270-280, 
1997. 


[27]  A.  Tolstoy,  K.  Smith,  and  N.  Maltsev,  “The  SWAM’99  Workshop  -  an  Over¬ 
view,”  Jowma/  of  Computational  Acoustics,  Vol.  9,  No.  1,  pp.  1-16,  2001. 

[28]  K.  B.  Smith,  “Computing  the  Influenee  of  Doppler  Due  to  Souree/Reeeiver  Mo¬ 
tion  in  a  Parabolie  Equation  Model,”  Journal  of  Computational  Acoustics,  Vol. 
10,No.  3,  pp.  295-309,  2002. 

[29]  G.  W.  Stewart,  “An  Updating  Algorithm  for  Subspaee  Traeking,”  IEEE  Transac¬ 
tions  on  Signal  Processing,  Vol.  40,  No.  6,  pp.  1535-1541,1992. 

[30]  G.  H.  Golub,  and  C.  F.  VanLoan,  Matrix  Computations,  3’^‘*  ed.,  John  Hopkins 
University  Press,  Baltimore,  1996. 


93 


[31]  G.  W.  Stewart,  Introduction  to  Matrix  Computations,  Academic  Press  Inc.,  Lon¬ 
don,  1973. 

[32]  Charles  W.  Therrien,  Discrete  Random  Signals  and  Statistical  Signal  Processing, 
Prentice-Hall,  Upper  Saddle  River,  New  Jersey,  1992. 

[33]  G.  W.  Stewart,  Matrix  Algorithms.  Volume  1:  Basic  Decompositions,  SIAM, 
Philadelphia,  1998. 


94 


INITIAL  DISTRIBUTION  LIST 


1 .  Defense  Technieal  Information  Center 
Ft.  Belvoir,  Virginia 

2.  Dudley  Knox  Library 
Naval  Postgraduate  Sehool 
Monterey,  California 

3.  Chairman,  Code  EC 

Eleetrieal  and  Computer  Engineering  Department 
Naval  Postgraduate  School 
Monterey,  California 

4.  Chairman,  Code  EA 

Engineering  Acoustics  Academic  Committee 
Naval  Postgraduate  School 
Monterey,  California 

5.  Prof  Roberto  Cristi,  Code  EC/Cx 

Electrical  and  Computer  Engineering  Department 
Naval  Postgraduate  School 
Monterey,  California 

6.  Mr.  Joe  Rice,  Code  PH/Rj 
Physics  Department 
Naval  Post  Graduate  School 
Monterey,  California 


95 


